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ABSTRACT 


The  primary  goal  of  this  work  is  to  use  mathematical  modelling  to  assist  in  defin¬ 
ing  the  operational  limits  of  the  Australian  Army  CH-47D  Chinook  when  carrying 
mixed  density  slung  loads.  This  report  presents  the  first  phase  in  the  program:  the 
development  of  a  simple  helicopter  slung-load  model  for  simulation  and  analysis  of  the 
system  dynamics. 

General  system  equations  of  motion  are  obtained  from  the  Newton-Euler  equations 
in  terms  of  generalized  coordinates  and  velocities.  The  system  is  partitioned  into 
coordinates  such  that  the  motion  due  to  cable  stretching  is  separated  from  that  due  to 
rigid-body,  coupled  dynamics.  In  the  formulation  used,  the  constraint  forces  appear 
explicitly  and  a  solution  to  the  resultant  generalized  accelerations  can  be  determined 
by  modelling  the  cable  as  a  simple  spring.  An  inelastic  solution  is  also  possible  by 
nulling  the  stretching  coordinates  to  obtain  a  relation  for  the  suspension  forces.  The 
system  equations  are  also  extended  for  the  multiple  load  case. 

The  model  is  verified  by  imposing  certain  constraints  in  order  to  approximate  a 
simple  pendulum  system  and  then  comparing  its  behaviour  with  analytical  results. 
Various  configurations  of  the  complete  helicopter  slung-load  system,  based  on  the 
CH-47B  Chinook  carrying  standard  military  containers,  are  then  examined  in  an  in¬ 
vestigation  of  the  open-loop  characteristics.  In  the  investigation,  several  parameters 
such  as  the  helicopter-load  mass  ratio,  suspension  configuration,  and  number  of  loads 
are  varied  and  the  resulting  system  modes  examined.  A  number  of  simulations  are 
also  presented  which  demonstrate  the  characteristic  behaviour  of  such  systems. 
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Mathematical  Modelling  of  Helicopter  Slung-Load  Systems 


EXECUTIVE  SUMMARY 

In  the  past,  the  operations  of  helicopters  carrying  externally  slung  loads  has  often 
been  limited  and,  in  some  cases,  seriously  hindered  by  stability  and  control  problems. 
Several  incidences  have  been  reported  by  the  Australian  Army  alone,  in  which  possible 
aerodynamic  excitation  or  dynamic  instability,  resulting  in  uncontrollable  oscillations,  has 
forced  premature  release  of  the  load.  Hence,  the  main  goal  of  the  study  is  to  develop  a 
comprehensive  helicopter  slung-load  model  which  will  provide  a  better  understanding  of 
the  system  dynamics  and  various  effects  involved.  Furthermore,  there  is  a  requirement  for 
the  carriage  of  multiple  loads  of  varying  type,  which  has  not  previously  been  investigated 
in  this  manner.  This  report  presents  the  first  phase  in  the  program:  the  development  of 
a  simple  helicopter  slung-load  model  for  simulation  and  analysis  of  the  system  dynamics. 
Following  this,  a  full  nonlinear  flight  dynamic  model  is  to  be  developed,  incorporating 
additional  detail,  such  as  the  automatic  flight  control  system,  rotor  wake  effects,  load 
aerodynamics,  and  sling  elasticity. 

The  simulation  model  used  for  this  first  phase  of  work  is  based  on  the  helicopter  slung 
load  system  first  introduced  at  NASA  Ames  Research  Center.  In  this  formulation,  the  gen¬ 
eral  system  equations  of  motion  are  obtained  from  the  Newton-Euler  equations  in  terms  of 
generalized  coordinates  and  velocities.  Using  the  explicit  constraint  method,  the  system  is 
then  partitioned  into  coordinates  such  that  the  motion  due  to  cable  stretching  is  separated 
from  that  due  to  rigid-body,  coupled  dynamics.  As  a  consequence,  the  constraint  forces 
appear  explicitly  and  a  solution  to  the  resultant  generalized  accelerations  is  determined  by 
assuming  a  simple  spring  model  for  the  cable.  It  is  also  possible  to  obtain  a  solution  to  the 
inelastic  approximation  by  nulling  the  stretching  coordinates,  yielding  an  explicit  relation 
for  the  suspension  forces.  The  result  is  computationally  more  efficient  than  the  conven¬ 
tional  formulation  and  is  readily  integrated  with  the  elastic  suspension  model.  Another 
benefit  of  the  formulation  is  that  it  is  easily  applied  to  complex  multiple  body  systems, 
and  in  the  current  work,  the  system  equations  are  extended  for  the  case  of  multiple  loads. 
All  code  development  has  been  done  in  the  MATLAB  numerical  computing  environment, 
which  provides  a  high-performance  language  amenable  to  modelling  and  simulation  type 
work.  The  main  functions  used  in  the  simulation  and  analysis  have  been  included  in  the 
report. 

The  model  is  verified  by  imposing  certain  constraints  in  order  to  approximate  a  simple 
pendulum  system  and  then  comparing  its  behaviour  against  both  analytical  solutions 
and  previously  documented  numerical  results.  A  complete  helicopter  slung-load  system, 
based  on  a  CH-47B  Chinook  carrying  a  standard  military  container,  is  then  examined  in 
an  investigation  of  the  open  loop  characteristics.  In  this  simple  model,  neither  the  load 
aerodynamics  nor  the  effect  of  rotor  downwash  on  the  load  is  taken  into  account.  Several 
parameters  such  as  the  helicopter-load  mass  ratio,  suspension  configuration,  and  number  of 
loads  are  varied  and  the  resulting  system  modes  examined.  In  order  to  extract  these  modes, 
a  linearised  form  of  the  model  is  first  obtained  by  numerical  approximation  of  the  partial 
derivatives.  For  the  configurations  examined,  an  increase  in  the  load  mass  was  generally 
found  to  have  a  mild  destabilising  effect,  particularly  in  the  lateral  axes.  A  number  of 
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simulations  were  also  run  in  order  to  demonstrate  the  dynamic  characteristics  of  several 
different  configurations.  From  the  response  data,  the  oscillatory  modes  in  longitudinal 
and  lateral  axes  were  identified,  including  the  unstable  phugoid  and  pendulum  modes. 
In  all  of  the  simulations  examined,  the  cable  tension  was  found  to  reach  a  maximum  of 
approximately  1.5  times  the  static  load. 
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A 
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c 

d 
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fg 
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9 
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m 
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n 
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eg  velocities  and  angular  velocities  of  the  n  bodies 

position  and  velocity  vectors  relative  to  inertial  space;  appended  num¬ 
bers  indicate  specific  locations  or  line  segments  joining  points;  the  su¬ 
perscript  star  denotes  the  eg  of  the  body 

suspension  force  parameters 

vector  of  cable  tensions  in  cable  Cj  ;  j  —  1, 2, . . .  ,m 

transformation  of  physical  vectors  from  frame  Fa  to  Fb  all  transfor¬ 
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generalised  velocity  coordinates  for  the  unconstrained  system 
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generalised  velocity  coordinates  defining  the  inelastic  and  elastic  com¬ 
ponents  of  the  system  motion,  respectively 

body  axes  configuration  velocity 

transformation  matrix  between  angle  rates  and  angular  velocities 
od  for  body  Bj 

vector  of  kinematic  accelerations  from  Euler’s  equations  for  each  body 

rigid  body  Euler  angle  triplet  and  inertial  angular  velocity;  appended 

numbers  indicate  specific  body 

vector  of  control  inputs  for  each  body 

moment  per  unit  tension  of  cable  Cj  on  body  Bj 

x,  y,  and  z  velocity  components  in  body  axes 
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aerodynamic  forces  along  each  corresponding  body  axes 

stability  and  control  derivatives  in  X  formulated  with  respect  to  the 
variables  defined  above 

stability  and  control  derivatives  in  Y 
stability  and  control  derivatives  in  Z 

aerodynamic  moments  about  each  corresponding  body  axes 

stability  and  control  derivatives  in  L 

stability  and  control  derivatives  in  M 

stability  and  control  derivatives  in  N 

moments  of  inertia  about  x,  y,  and  z  axes 

product  of  inertia  in  x  —  z  plane 

dynamic  pressure 


Operators 
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0. 

()r 

0* 

S(Va) 


dot  and  cross  product  operators  for  physical  vectors 
physical  vector  given  by  its  coordinates  in  the  frame  Ta 
transpose  of  () 

quantity  associated  with  eg  of  a  rigid  body  in  the  system 
skew-symmetric  matrix  representing  cross-product  operation  for  vec¬ 
tors,  V  referred  to  Ta 
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1  Introduction 


The  operations  of  helicopters  carrying  externally  slung  loads  has  often  been  limited 
and,  in  some  cases,  seriously  hindered  by  stability  and  control  problems.  Several  incidences 
have  been  reported  by  the  Australian  Army  alone  in  which  aerodynamic  excitation  or 
dynamic  instability,  resulting  in  uncontrollable  oscillations,  has  forced  premature  release 
of  the  load. 

A  program  was  consequently  initiated  within  the  Defence  Science  and  Technology 
Organisation  (DSTO)  to  use  computer  modelling  and  simulation  to  assist  in  defining  the 
operational  limits  of  the  Australian  Army  CH-47D  Chinook  when  carrying  slung  loads. 
The  first  phase  in  this  program  has  entailed  the  development  of  a  simple  helicopter  slung- 
load  model  for  simulation  and  analysis  in  order  to  provide  a  better  understanding  of  the 
system  dynamics  and  various  effects  involved.  The  results  from  this  work  are  presented 
here.  In  the  second  phase,  a  more  comprehensive  helicopter  model  is  to  be  integrated 
into  the  multi-body  system.  The  simulation  model  will  incorporate  additional  detail,  such 
as  the  automatic  flight  control  system,  load  aerodynamics,  rotor  wake  effects,  and  sling 
elasticity.  Furthermore,  there  is  a  requirement  to  model  the  dynamics  of  the  helicopter 
with  multiple  slung  loads  of  varying  mass  and  aerodynamic  properties,  which  has  not 
previously  been  investigated  in  this  manner. 

In  this  report,  Section  2  presents  a  broad  overview  of  prior  research  into  helicopter 
slung-load  systems.  The  three  areas  of  analytical  studies,  experimental  testing  and,  flight 
simulation  are  covered. 

Section  3  introduces  the  equations  of  motion  for  a  generalised  helicopter  slung-load 
model.  Both  inelastic  and  elastic  formulations  axe  included.  The  overall  simulation  pro¬ 
cedure  is  then  explained  in  general  terms.  The  full  set  of  equations  which  constitute  the 
simulation  model  are  listed  in  Appendix  A. 

In  Section  4,  the  results  obtained  from  the  analysis  of  a  simple  pendulum-type  model 
are  presented.  These  include  a  comparison  of  the  characteristic  system  modes  and  longi¬ 
tudinal  and  lateral  simulations  of  two  different  configurations. 

In  Section  5,  the  dynamics  of  the  helicopter  slung-load  model  are  discussed.  The  modes 
of  the  linearised  system  are  displayed  over  a  range  of  mass  and  sling  configurations  and 
the  characteristic  behaviour  of  a  Super  Stallion  CH-53D  cargo  helicopter  with  MILVAN 
container  load  are  compared  against  previous  results. 

Finally,  some  concluding  remarks  are  drawn  and  proposals  for  further  research  made. 
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2  Background 

There  has  been  a  small  but  significant  amount  of  work  done  in  investigating  the  be¬ 
haviour  and  control  of  helicopter  slung-load  systems.  This  work  can  be  roughly  divided 
into  three  veins:  analytical  studies,  experimental  testing,  and  flight  simulation  and  au¬ 
tomatic  control.  In  the  early  years,  most  of  this  effort  was  concentrated  in  analytical 
studies,  including  various  control  designs.  Experimental  testing  in  flight  and  wind  tunnel 
was  mainly  limited  to  the  establishment  of  operational  limits  based  on  gross  aerodynamic 
instabilities.  With  the  advancement  of  digital  computing,  however,  came  the  ability  to 
perform  dynamic,  piloted  simulations  in  real  time  and  develop  more  complex  models. 
This  increase  in  the  complexity  of  the  system  led  to  a  requirement  for  better  aerodynamic 
models  —  for  both  helicopter  and  load  —  and  the  emphasis  in  experimental  work  shifted 
accordingly.  These  days,  there  is  still  some  analytical  work  being  conducted,  particularly 
in  the  design  of  automatic  control  systems,  but  much  of  the  research  is  now  undertaken 
using  simulation. 


2.1  Analytical  Studies 

Some  of  the  earliest  studies  in  helicopter  slung-load  behaviour  was  carried  out  by 
Lucassen  and  Sterk  [1]  in  1965  to  provide  a  better  understanding  of  the  dynamics  and 
indicate  means  to  avoid  undesirable  effects,  since  they  were  known  to  cause  a  reduction 
in  the  helicopter’s  operational  capabilities.  In  this  work,  the  equations  of  motion  were 
restricted  to  the  vertical  plane  with  three  degrees  of  freedom  (DOF),  and  load  aerody¬ 
namics  was  not  included.  Later,  Dukes  [2]  used  a  similar  approximation  to  examine  the 
modes  of  the  system  in  the  frequency  domain  and  explore  various  feedback  and  open-loop 
control  systems  for  damping  the  pendulous  helicopter-load  motion.  Cliff  and  Bailey  [3] 
also  used  a  simplified  model  of  a  helicopter  with  a  singly  tethered  load,  which  neglected 
all  aerodynamics  other  than  drag.  In  their  formulation,  the  nonlinear  equations  of  motion 
were  linearised  about  a  steady  level  flight  condition,  and  then  the  resulting  perturba¬ 
tion  equations  were  separated  into  longitudinal  and  lateral- directional  sets.  Although  the 
authors  were  able  to  make  some  inferences  regarding  the  stability  effects  of  various  param¬ 
eters,  such  as  mass  ratio  and  tether  length,  they  suggested  that  more  complete  dynamical 
models  needed  investigation.  An  attempt  to  increase  the  fidelity  of  the  aerodynamic  load 
model  was  made  by  Feaster  [4]  and  Feaster  et  al  [5]  using  an  experimentally  determined 
yaw-damping  coefficient  in  a  linearised  small  perturbation  stability  analysis,  which  con¬ 
sidered  both  single-cable  and  two-cable  tandem  suspension  systems.  The  results  agreed 
well  with  the  previous  model  and  full-scale  tests  and  also  demonstrated  that  the  two-cable 
tandem  suspension  system  offered  a  satisfactory  means  of  transporting  the  standard  cargo 
container  examined.  Around  the  same  time,  Prabhakar  and  Sheldon  [6]  and  Prabhakar 
[7]  undertook  a  theoretical  study  of  a  Westland  Sea  King  helicopter  carrying  a  standard 
cargo  container  on  a  two  point  longitudinal  suspension.  Again,  the  aerodynamic  stability 
derivatives  used  in  the  model  were  determined  through  experiment  and  it  was  found  that 
the  pitch  and  yaw  rate  derivatives  were  strongly  destabilising. 

Over  the  following  several  years,  Nagabhushan  [8,  9]  and  Nagabhushan  and  Cliff  [10] 
produced  several  reports  on  the  dynamics,  stability,  and  control  of  helicopters  carrying 
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externally  suspended  loads.  Several  mathematical  models  of  varying  order  were  developed 
to  describe  the  dynamics  of  such  systems.  However,  unlike  much  of  the  previous  analytical 
work,  which  was  based  on  the  Newton-Euler  equations,  the  nonlinear  formulation  was 
derived  from  Lagrange’s  equations  for  general  dynamical  systems.  In  these  reports,  the 
low-speed  stability  characteristics  of  a  conventional  helicopter  with  an  external  sling  load 
on  a  single-point  suspension  were  investigated.  Typically,  towing  cable  length,  towed 
body-to- vehicle  mass  ratio,  and  load  factor  in  a  turn  were  found  to  affect  the  stability  of 
the  aircraft  and  its  sling  load. 

In  1986,  Ronen  [11]  and  Ronen  et  al  [12]  developed  a  new  model  for  a  helicopter 
carrying  a  sling  load  on  a  single  point  suspension  in  order  to  improve  on  the  existing 
dynamic  models  and  investigate  the  open  loop  characteristics  of  the  system.  For  the  first 
time,  the  model  took  into  account  the  effects  of  rotor  downwash  on  the  load  and  the 
unsteady  aerodynamics  of  bluff-body  type  loads.  The  nonlinear  equations  of  motion  were 
derived  and  then  separated  into  two  sets:  the  nonlinear  trim  equations  and  the  linearised 
equations  for  small  perturbation  about  the  equilibrium.  More  recently,  Curtiss  [13]  derived 
a  full  set  of  equations  for  the  twin  lift  system,  linearised  about  a  hover  trim  condition. 


2.2  Experimental  Testing 

Although  there  has  been  a  substantial  amount  of  experimental  work  in  determining 
the  aerodynamic  behaviour  of  various  slung  loads  using  wind-tunnels,  little  has  been  done 
through  full  scale  flight-testing.  In  1968,  Gabel  and  Wilson  [14]  presented  the  results  of 
an  extensive  program  of  simulation,  wind-tunnel,  and  flight  tests,  which  were  conducted 
to  assist  in  solving  the  problems  of  sling  load  vertical  bounce,  sling-leg  web  flapping,  and 
aerodynamic  yaw  instability.  Some  years  later,  Hone  [15]  utilised  the  data  from  actual 
flight  tests  on  a  Sikorsky  CH-54  heavy-lift  helicopter  to  investigate  the  validity  of  a  model 
developed  by  Briczinski  and  Karas  [16].  The  aim  of  this  work  was  to  explore  phenomena 
associated  with  the  carriage  of  externally  suspended  loads  on  helicopters,  and  to  establish 
more  reliable  strength  requirement  data  for  the  load  slings  and  their  interfaces. 

Other  work  in  experimental  testing  includes  those  presented  by  Kesler  et  al  [17],  Feaster 
[4],  Feaster  et  al  [5],  and  Matheson  [18]. 

2.3  Flight  Simulation  and  Control 

One  of  the  first  investigations  into  automatic  control  for  helicopters  with  slung  loads 
was  conducted  by  Wolkovitch  and  Johnston  [19]  in  1965.  The  single-cable  dynamic  model 
was  developed  in  a  straightforward  application  of  the  Lagrange  equations.  Abzug  [20] 
later  expanded  on  this  model  to  consider  the  case  of  two  tandem  cables.  However,  his 
formulation  was  based  on  the  Newton-Euler  equations  of  motion  for  small  perturbations, 
separated  into  longitudinal  and  lateral  sets.  Aerodynamic  forces  on  the  cables  and  the 
load  were  neglected,  as  were  the  rotor  dynamic  modes. 

In  recognition  of  suspension- related  problems  encountered  with  the  carriage  of  external 
cargo  by  helicopters,  the  US  Army  in  1970  initiated  a  program  aimed  at  the  establishment 
of  design  criteria  for  sling  members  and  hard-points.  This  program,  as  well  as  many 
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subsequent  investigations,  were  undertaken  by  the  Eustace  Directorate,  US  Army  Air 
Mobility  Research  and  Development  Laboratory  (USAAMRDL).  Part  of  the  first  phase  in 
the  contract,  reported  by  Briczinski  and  Karas  [16],  involved  the  computerised  simulation 
of  a  helicopter  and  external  load  in  real  time  with  a  pilot  in  the  loop.  Load  aerodynamics 
were  incorporated  into  the  model,  as  well  as  rotor-downwash  effects  in  hover.  Soon  after 
this  program,  Liu  [21]  conducted  an  extensive  study  to  select  the  best  technical  approaches 
for  stabilising  a  wide  spectrum  of  externally  slung  helicopter  loads  at  forward  speeds. 
The  simulation  model  used  extended  that  of  Abzug  [20]  to  include  load  aerodynamics. 
Several  stabilisation  systems  were  evaluated  using  a  moving-base  simulator  and  of  those, 
an  electronic  system  providing  rate  and  acceleration  inputs  to  the  helicopter’s  stability 
augmentation  system  (SAS)  was  favoured.  The  design  and  assessment  of  automatic  control 
systems  continued  with  Asseo  and  Whitbeck  [22]  in  their  paper  on  the  control  requirements 
for  sling-load  stabilisation.  Linearised  equations  of  motion  of  the  helicopter,  winch,  cable, 
and  load  complex  were  developed  for  a  variable  suspension  geometry  and  were  then  used 
in  conjunction  with  modern  control  theory  to  design  several  control  systems  for  each  type 
of  suspension.  The  next  year,  Gera  and  Farmer  [23]  examined  the  feasibility  of  stabilising 
external  loads  by  means  of  controllable  fins  attached  to  the  cargo.  In  their  simple  linear 
model  representing  the  yawing  and  the  pendulous  oscillations  of  the  slung-load  system,  it 
was  assumed  that  the  helicopter  motion  was  unaffected  by  the  load. 

Following  the  first  program  of  work  sponsored  by  the  Eustace  Directorate,  a  further 
study  to  define  important  flight  control  system  design  and  handling  qualities  criteria  for 
moving  loads  slung  beneath  tandem-rotor  helicopters  was  conducted  by  Kesler  et  al  [17]. 
It  included  theoretical  analyses,  acquisition,  and  evaluation  of  both  wind  tunnel  and  flight 
test  data,  analysis  of  various  problems,  and  the  actual  flight  simulation  of  a  Model-347 
advanced  tandem-rotor  helicopter  with  an  external  load.  Another  program  under  the  same 
sponsorship,  investigated  by  Alansky  et  al  [24],  looked  into  the  quantitative  limitations 
of  the  CH-47  helicopter  performing  terrain  flying  with  external  loads.  The  simulation 
used  in  this  investigation  comprised  a  fully  coupled  total  force  and  moment  model  and 
an  alternative  method  of  load-control  named  the  Active  Arm  External  Load  Stabilisation 
System  (AAELSS). 

Some  time  later,  a  generalised  real  time,  piloted,  visual  simulation  of  a  single  rotor 
helicopter,  suspension  system,  and  external  load  was  developed  by  Shaughnessy  et  al  [25], 
and  subsequently  validated  for  the  full  flight  envelope  of  a  CH-54  helicopter  and  cargo 
container.  The  mathematical  model  described  used  modified  nonlinear  classical  rotor 
theory  for  both  the  main  rotor  and  tail  rotor,  nonlinear  fuselage  aerodynamics,  an  elastic 
suspension  system,  nonlinear  load  aerodynamics,  and  a  load-ground  contact  model. 

In  1980,  Sampath  [26]  completed  his  PhD  dissertation  on  the  dynamics  of  a  tandem- 
rotor  helicopter  slung-load  system,  which  involved  modelling  and  simulation  as  well  as 
experimental  wind-tunnel  tasks.  In  his  formulation,  Lagrange’s  equations  were  used  to 
write  the  equations  of  motion  and  were  divided  into  two  sets:  one  for  the  towing  vehicle 
and  the  other  for  the  slung  load.  The  cables  of  the  sling  were  modelled  as  massless  linear 
springs  with  viscous  damping  and  no  aerodynamic  properties.  The  aerodynamic  models 
for  the  helicopter  and  load  were  both  implemented  using  tabulated  static  data.  Some 
years  later,  a  full  nonlinear  simulation  model  of  the  CH-47B  helicopter,  developed  by  the 
Boeing  Vertol  Company,  was  adapted  for  use  in  the  NASA  Ames  Research  Center  (ARC) 
simulation  facility  by  Weber  et  al  [27].  The  mathematical  model  developed  was  based 
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on  a  total  force  approach  in  six  rigid-body  DOF  along  with  the  option  for  an  externally 
suspended  load  in  three  DOF.  The  aerodynamic  models  were  also  quite  comprehensive, 
including  steady-state  rotor  flapping  and  load  aerodynamic  effects. 

More  recently,  research  into  the  automatic  control  of  load  dynamics  has  continued 
with  Raz  et  al  [28]  in  an  investigation  of  an  active  aerodynamic  Load  Stabilisation  System 
(LSS)  for  a  helicopter  sling-load  system.  The  general  theoretical  model  used  was  based  on 
previous  work  [11]  for  configurations  with  a  single  suspension  point. 

Some  of  the  most  recent  work  in  the  simulation  of  helicopter  slung-load  systems  has 
been  conducted  by  Cicolani  and  Kanning  [29,  30]  and  Cicolani  et  al  [31]  at  the  NASA 
ARC.  In  this  work,  the  general  simulation  equations  were  derived  for  the  motion  of  slung- 
load  systems  consisting  of  several  rigid  bodies  connected  by  straight-line  cables  or  links, 
assumed  to  be  either  elastic  or  inelastic.  A  formulation  for  the  general  system  was  obtained 
from  the  Newton-Euler  rigid-body  equations  with  the  introduction  of  generalised  velocity 
coordinates.  The  same  approach  for  simulating  helicopter  slung-load  dynamics  has  been 
adopted  in  the  current  work  at  AMRL.  In  addition,  the  equations  of  motion  have  been 
extended  to  the  case  of  multiple  loads  with  disparate  mass  and  aerodynamic  characteristics 
and  sling  configurations. 


2.4  Surveys  and  Overviews 

In  his  book  on  the  dynamics  of  helicopter  flight,  Saunders  [32]  devoted  one  section 
to  piloting  problems  associated  with  the  carriage  of  external  loads.  The  normal  trim 
state  and  origins  of  load  oscillation  were  first  examined  using  a  simple  helicopter-load 
mass  model.  Various  problems,  including  uncontrollable  load  oscillations,  aerodynamic 
excitation,  active  control,  and  poor  visibility  were  then  discussed  in  a  broad  context. 

In  1976,  Matheson  [18]  compiled  a  review  of  the  developments  and  data  concerning 
the  operation  of  helicopters  with  slung  loads.  The  report  focused  on  the  problems  of 
aerodynamic  instability,  vertical  bounce,  and  sling-leg  flapping.  Methods  for  reducing 
these  instabilities  and  procedures  for  extending  the  operating  limits  of  a  helicopter  with 
different  types  of  slung  loads  were  also  discussed. 

Around  the  same  time,  Shaughnessy  and  Pardue  [33]  conducted  a  survey  of  the  heli¬ 
copter  sling  load  accident/incident  records  provided  by  various  US  organisations  for  the 
period  from  1968  to  1974.  From  the  data,  the  highest  percentage  of  accidents  occurred 
during  hover,  and  it  was  therefore  concluded  that  hovering  was  the  most  critical  sling 
load  flight  operation.  Furthermore,  the  accidents  and  incidents  caused  by  swinging  loads 
during  cruise  were  generally  much  less  severe  than  the  hover  mishaps. 
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3  System  Representation 

3.1  Description  of  System 

Helicopter  slung-load  systems  fall  into  a  class  of  multibody  dynamic  systems  consisting 
of  two  or  more  rigid  bodies  connected  by  massless  links.  The  links  can  be  considered  either 
elastic  or  inelastic,  although  the  rigid-body  assumption  excludes  any  helicopter  or  load 
elastic  modes.  Typically,  the  system  is  characterised  by  the  configuration  geometry,  mass, 
inertia,  and  aerodynamic  behaviour  of  both  helicopter  and  load,  as  well  as  the  elastic 
properties  of  the  links. 

In  general  terms,  the  system  of  interest  consists  of  a  single  helicopter  supporting  one  or 
more  loads  by  means  of  some  suspension.  Several  examples  of  the  various  configurations 
under  consideration  are  illustrated  in  Figure  1.  The  model  is  comprised  of  n  rigid  bodies, 
Bl,B2,...,Bn,  with  m  straight-line  links  supporting  a  single  force  in  the  direction  of  the 
link.  For  cables,  this  is  strictly  a  tensile  force  —  cable  collapse  is  not  considered.  If  the 
links  are  modelled  as  inelastic,  c  (<  m)  holonomic  constraints1  are  imposed  on  the  motion 
of  the  bodies  and  the  system  has  d  =  6n  —  c  DOF.  If  the  links  are  modelled  as  elastic, 
there  are  6 n  DOF. 

In  the  model  used,  a  number  of  simplifying  assumptions  were  made.  These  included 
the  exclusion  of  cable  aerodynamics  and  rotor-downwash  effects.  Furthermore,  load  aero¬ 
dynamics  have  been  neglected  for  this  initial  stage  of  the  work.  Despite  these  limitations, 
the  system  defined  above  has  proven  adequate  for  simulation  studies  [31]  in  which  the 
low-frequency  behaviour  is  of  primary  interest. 


3.2  Generalized  Equations  of  Motion 

The  simulation  model  used  for  this  first  stage  of  work  was  based  on  the  helicopter 
slung-load  system  introduced  by  Cicolani  and  Kanning  [29].  In  this  formulation,  the 
general  system  equations  of  motion  are  obtained  from  the  Newton-Euler  equations  in 
terms  of  generalised  coordinates  and  velocities.  Following  the  explicit  constraint  method, 
which  utilises  d’Alembert’s  principle,  the  system  is  partitioned  into  coordinates  such  that 
the  motion  due  to  cable  stretching  is  separated  from  that  due  to  rigid-body,  coupled 
dynamics.  As  a  consequence,  the  constraint  forces  appear  explicitly  and  a  solution  to  the 
resultant  generalised  accelerations  is  determined  by  assuming  a  simple  spring  model  for 
the  cable. 

It  is  also  possible  to  obtain  a  solution  to  the  inelastic  approximation  by  nulling  the 
stretching  coordinates  to  obtain  an  explicit  relation  for  the  suspension  forces.  The  result 
is  computationally  more  efficient  than  conventional  procedures  and  is  readily  integrated 
with  the  formulation  for  elastic  suspension.  Another  benefit  of  the  formulation  is  that 
it  is  easily  applied  to  complex,  multiple  body  systems,  as  in  the  current  work.  To  date, 
all  code  development  has  been  done  in  the  MATLAB  numerical  computing  environment, 

holonomic  constraints  represent  excess  coordinates  which  are  independent  and  can  be  eliminated 
through  equations  of  constraint. 
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which  provides  a  high-performance  language,  amenable  to  modelling  and  simulation  type 
work. 


Figure  1:  Single  Point  Slung-Load  Configurations 


The  Newton-Euler  equations  of  motion  for  a  system  of  n  rigid  bodies  can  be  expressed 
in  six  DOF  as 


7m'<7iv  -I-  FAix  +  FCi jv  —  miVi*N  =  0 

MAii  +  MCi{  —  Jiujii  —  S(uii)Jiuii  =  0 


*  =  1,2,. 


In  this  expression,  the  first  set  of  equations  represents  the  balance  of  translational 
forces,  where  the  subscript  N  denotes  the  inertial  axes.  The  second  set  represents  the  sum 
of  moments  about  each  body’s  eg,  where  the  subscript  i  denotes  the  corresponding  body 
axes.  Both  equations  are  comprised  of  several  terms  including  the  forces  and  moments  due 
to  gravity,  aerodynamics,  and  inertia.  The  first  term,  mig^,  is  the  gravity  force  acting 
through  each  eg,  FAi jv  and  MAiz  are  the  aerodynamic  forces  and  moments  respectively, 
and  FCifj  and  MCil  the  cable  forces  and  moments  respectively.  The  terms  miVi*N  and 
Ji  uii  constitute  the  inertial  reaction  of  each  body,  and  S (c oii)  J i  toil  is  the  moment  induced 
by  the  Coriolis  effect. 

It  is  convenient  to  write  these  equations  as  a  single  expression  in  matrix  form.  Denoting 
the  configuration  position  vector  as  r  and  the  configuration  velocity  as  v  for  the  system, 


where  the  rigid-body  eg  positions  Ri*N  =  [  xt  yi  zx  ] T  and  Euler  angles  ai  =  [fa  Qi  (pi  ]  - 
The  corresponding  velocities  Vi*N  =  Ri*N  =  [  x%  yx  i;  ]  and  the  angular  rates  uii  = 

[  Pi  Qi  ri  ]T  are  related  to  the  Euler  angle  rates  via  the  transformation 

uii  —  Wiidii  (3) 


where 


1  0  —  sin  6i 

Wii  =  0  cos  fa  sin  fa  cos  6X 

0  -  sin  fa  cos  fa  cos  9{ 


(4) 
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Using  fg,  fa,  fc,  and  /*,  for  the  combined  force- moment  vectors  due  to  gravity, 
aerodynamics,  cable  suspension,  and  inertia  plus  Coriolis  effects,  the  equations  of  motion 
can  be  written 


fg  +  fa  +  fc  +  f*  =  0 


(5) 


where 


and 

where 


'  ml  gN  ' 

'  fain  ' 

'  FC1N  ‘ 
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0 

J1 
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Jn 

S(unn)  Jnu>nn 

(6) 


(7) 


(8) 


Here,  D  is  a  block-diagonal  matrix  comprising  masses  and  inertias  along  the  main 
diagonal,  v  is  the  configuration  acceleration  and  the  matrix  X  contains  the  Coriolis  terms 
due  to  the  use  of  rotational  coordinates  in  body-axes 

In  order  to  derive  a  set  of  simulation  equations  for  the  system,  a  solution  to  the  equa¬ 
tions  of  motion  described  above  must  be  found;  that  is,  an  expression  for  the  configuration 
acceleration  in  terms  of  the  system  states  and  applied  forces.  For  the  helicopter  slung-load 
system  under  consideration,  it  is  useful  to  first  formulate  a  set  of  generalized  coordinates 
and  velocities  which  describe  the  motion  of  the  inelastic  system  and  the  effect  of  cable 
stretching  as  two  distinct  subsets.  An  inelastic  system  with  n  bodies  and  c  constraints  will 
have  d  =  6n  -  c  DOF.  The  cable  constraints  on  the  helicopter  slung-load  system  can  be 
considered  holonomic.  In  addition,  for  the  following  system,  the  constraints  will  be  posed 
as  time  invariant.  The  special  cases  of  cable  winching  and  attachment-point  movement 
are  not  considered. 


The  system  can  be  partitioned  according  to  the  6 n  generalized  position  coordinates  as 
follows: 

’  gi 
A 


<7  = 


(9) 


where  ql  is  the  list  of  d  position  coordinates  for  a  system  with  inejastic  suspension  and 
A  axe  the  c  coordinates  which  describe  the  variation  in  cable  length  due  to  stretching. 
The  configuration  velocity  can  be  expressed  as  a  linear  function  of  the  generalized  velocity 
coordinates 

v  =  Au  (10) 
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where 


u 


ul 

A 


A 


A1 

L 


Differentiating  Equation  (10)  and  substituting  v  into  Equation  (7)  yields 

f*  =  -DAu  -  DAu  -  X 


(11) 


(12) 


Now,  replacing  /*  in  Equation  (5),  the  following  simplified  version  of  the  equations  of 
motion  can  be  written: 

fo  +  fc-  DAu  =  0  (13) 

where  the  vector  fo  is  the  sum  of  all  external  forces  and  inertial  coupling  terms,  i.e. 

fo  =  fg  +  fa  -  DAu  -  X  (14) 

Since  the  system  has  been  specified  in  terms  of  its  generalized  coordinates,  A  is  a  square 
6n  x  6n  nonsingular  matrix  and  a  solution  for  the  generalized  acceleration  coordinates 
exists.  Prom  Equation  (13),  the  acceleration  equation  is 


ii  =  A  'D  '[fo  +  fc] 


(15) 


It  should  be  noted  here  that  the  inverse  matrix  A-1  simply  represents  the  relation 
u(v)  that  can  be  derived  analytically  from  the  kinematics.  Therefore,  it  is  unnecessary 
to  perform  a  costly  numerical  inversion  to  obtain  u.  In  partitioned  form,  the  acceleration 
equation  is 


ul 

'  AI1t  ' 

A 

AT 

D  '[fo  +  fc] 


(16) 


where  AI1T  and  AT  are  the  6 n  —  c  and  c  rows  of  A-1  which  define  the  relations  ul  and 
A  respectively.  Prom  the  first  set  of  equations  representing  the  inelastic  component,  the 
solution  for  ul  is 


ill  =  AI1T D  '[fo  +  fc] 


(17) 


The  second  set  of  equations  represents  the  elastic  component 

A  =  A  TD~'[fo  +  fc]  (18) 


An  alternate  formulation  for  the  accelerations  ul  can  be  obtained  by  first  differentiating 
the  expanded  form  of  the  configuration  velocity  from  Equation  (10)  as  follows: 

v  =  Alul  +  LA  (19) 

and  then  substituting  into  Equation  (7)  to  give 

/*  =  -DAlul  -  DAliil  -  D{LA  +  LA)  -  X  (20) 

If  fo  is  redefined: 

fo  =  fg  +  fa-  DAlul-  X  (21) 


9 


DSTO-TR-1257 


then  replacing  f*  in  Equation  (5)  and  solving  for  ul  produces 

ill  =  [AlTDAl]-lAlT[fo  -  D(LX  +  LX)  +  fc]  (22) 

The  last  step  required  in  determining  a  solution  for  the  generalized  accelerations  is  to 
calculate  the  constraint  force  fc.  For  a  system  with  c  constraints,  the  constraint  force  can 
be  expressed  as 

fc  =  Hs  (23) 

where  the  columns  of  the  matrix  H  are  configuration  vectors  and  rank  H  —  c.  The 
elements  of  the  vector  s  are  arbitrary  scalars  to  be  determined.  The  exact  form  of  this 
equation  and  its  solution  depend  on  whether  the  cables  are  considered  elastic  or  not. 


Elastic  System 

For  a  general  elastic  system  with  M  cable-body  attachments,  the  suspension  forces  on 
each  body  can  be  given  as  the  sum  of  forces  and  moments  applied  at  each  attachment 
point  on  that  body,  i.e. 

M 

fc  =  J2  hj  Tj  (24) 

3= 1 

In  this  formulation,  the  configuration  vector  hj  defines  the  force  and  eg  moment  due  to 
a  unit  load  at  the  jth  cable  attachment,  and  Tj  is  the  cable  tension.  Referring  to  Figure  2, 
kcj  and  (Ri*j  x  kcj)  denote  the  constraint  force  and  moment  per  unit  tension  respectively, 
on  body  Bi(j)  due  to  the  jth  cable  attachment.  The  vector  kcj  is  the  cable  direction 
outward  from  the  body,  and  Ri*j  =  (Rj  —  Ri*)  is  the  moment-arm  of  the  attachment 
point  about  the  eg. 


Figure  2:  Suspension  Forces  on  a  Rigid  Body 

The  tension  in  each  cable  is  given  by  a  simple  spring  model  as 

Tj  =  max{0,  Kj (tj  -  loj }  +  cjtj}  ;  j  =  1,2,...  ,M  (25) 

where  loj  and  tj  are  the  unloaded  and  instantaneous  cable  lengths,  and  Kj  and  Cj  are 
the  cable  spring  and  damping  constants. 
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Inelastic  System 

It  can  be  shown  that  the  columns  of  H  and  A  both  form  bases  of  the  same  linear  vector 
space  and  therefore  A  can  be  used  to  define  the  constraint  force,  i.e. 

fc  =  As  (26) 

where  the  vector  s  will  have  units  of  force  if  the  coordinates  A  are  lengths.  To  find  a 
solution  for  the  inelastic  system,  the  constraint  acceleration  is  set  to  zero,  i.e.  A  =  0. 
Substituting  into  Equation  (18)  gives 

0  =  ATD~l[fo  +  As]  (27) 

the  solution  to  which  is 

s  =  —[ATD~1A]~1ATD~1fo  (28) 


3.3  Simulation  of  Helicopter  Slung  Load  System 

Prior  to  executing  the  simulation,  several  components  must  be  customised  to  the 
helicopter  slung- load  configuration  of  interest. 

The  first  step  in  setting  up  the  simulation  equations  involves  determining  the  con¬ 
straints  of  the  inelastic  system  and  then  defining  the  generalised  velocity  coordinates  (ul, 
A).  Using  these  coordinates  in  kinematic  relations  for  the  system,  it  is  then  possible  to 
obtain  expressions  for  the  system  matrices  A,  A -1,  and  A.  The  selection  of  appropriate 
coordinates  is  case  specific;  however,  it  is  possible  to  choose  them  so  that  they  consist 
largely  of  natural  vectors.  In  most  applications,  including  that  discussed  in  this  report, 
u  is  comprised  of  the  eg  velocity  of  a  reference  body  (typically  the  helicopter),  the  cable 
velocities,  and  the  angular  velocities  of  all  bodies  including  both  helicopter  and  loads. 

Next,  an  appropriate  representation  for  the  suspension  cables  must  be  chosen.  For 
inelastic  cables,  fc  is  calculated  from  any  basis  of  the  constraint  force  space  A  and  the 
corresponding  constraint  force  parameters  s  as  in  Equation  (28). 

For  the  last  step,  the  aerodynamic  and  inertial  properties  of  both  the  helicopter  and 
loads  need  to  be  implemented  in  the  model.  For  most  rigid  bodies,  the  aerodynamic 
forces  and  moments  are  a  function  of  the  configuration  velocities  and  displacements  v 
and  r,  and  the  control  inputs  5.  Typically,  the  helicopter  aerodynamic  model  neglects 
position-dependent  and  acceleration-dependent  effects  such  as  interbody-ground  interfer¬ 
ence  and  unsteady  aerodynamics.  However,  these  are  often  secondary  in  nature  and  the 
resulting  model  is  adequate  for  simulation  under  most  conditions.  Aerodynamic  models 
for  loads,  which  are  generally  unsteady  and  of  much  higher  order,  axe  less  well  understood 
or  replicated. 

Once  the  system  has  been  configured,  the  dynamic  simulation  can  proceed.  First,  the 
initial  state  (u,r)  and  the  trim  state  ( uq,tq )  must  be  set.  Then  the  integration  loop  is 
started  and  the  following  steps  are  executed  in  sequence,  according  to  the  flow-diagram  of 
Figure  3: 
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Figure  3:  Simulation  Flow  Diagram 
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1.  Determine  the  aerodynamic  force  fa,  inertia,  and  Coriolis  forces  /*,  and  gravity 
force  fg,  all  in  inertial  axes.  Assuming  the  aerodynamic  model  is  written  in  body 
axes,  an  angular  transformation  will  be  required  for  fa.  The  configuration  velocity 
v  can  be  calculated  from  the  generalised  velocity  u  using  the  kinematic  matrix  A. 

2.  Sum  the  external  forces  fa,  /*,  and  fg  to  yield  the  configuration  force  fo. 

3.  Using  the  configuration  force  along  with  matrices  derived  from  the  current  state 
( u ,  r)  and  the  configuration  geometry,  solve  the  cable  force  fc  for  either  elastic  or 
inelastic  suspension  models. 

4.  Compute  the  generalised  accelerations  it  from  the  inverse  kinematic  matrix  A~l,  and 
the  inverse  mass  matrix  J9  . 

5.  Compute  the  velocity  r  from  the  configuration  velocity  and  inverse  transformation 
matrix  W  . 

6.  Apply  an  integration  step  to  predict  the  new  state  ( u ,  r)  and  repeat  the  sequence. 

At  this  early  stage,  all  code  development  has  been  done  in  the  MATLAB  [34]  numerical 
computing  environment,  which  provides  a  high-performance  language  amenable  to  mod¬ 
elling  and  simulation  work.  It  is  important  to  stress  that  this  pilot  simulation  was  not 
intended  to  run  in  real  time,  but  rather  produce  the  appropriate  output  for  subsequent  re¬ 
play  and  analysis.  At  a  later  stage,  the  code  will  be  ported  to  a  platform-specific  compiled 
language  suitable  for  piloted,  real  time  simulation. 

The  Helicopter  Slung-Load  Simulation  (HSLSIM)  program,  consists  of  several  modules. 
These  including  the  main  script,  integration  function,  differential  equation  solution,  aero¬ 
dynamic  model,  and  various  output  and  replay  functions.  The  simulation  is  run  through 
the  main  script,  which  generates  the  control  inputs,  configures  the  helicopter-load  system 
properties  (geometric  and  inertial),  sets  the  initial  system  state,  and  then  executes  the 
integration  function.  The  integration  function  ODE45  is  problem  independent  and  based 
on  an  algorithm  which  combines  4th  and  5th  order  Runge-Kutta  formulas  for  ordinary  dif¬ 
ferential  equations.  It  requires  a  function  tailored  to  the  problem  at  hand,  which  provides 
a  point  solution  to  the  differential  equation.  For  the  helicopter  slung-load  simulation,  this 
function  represents  the  core  of  the  code  and  implements  much  of  the  above  flow  diagram. 
The  aerodynamic  models  for  both  helicopter  and  loads  axe  called  from  within  this  func¬ 
tion.  They  can  be  as  simple  or  as  complex  as  desired,  but  must  output  total  force  and 
moment  variables.  Hence,  if  small-perturbation  aerodynamic  models  axe  to  be  used,  they 
must  be  augmented  with  the  corresponding  trim  forces  and  moments.  The  cable  elastic 
model  can  also  be  implemented  as  a  separate  function,  although  this  was  not  done,  since 
the  spring-damper  model  is  fairly  standard  and  easily  included  in  the  solution  function. 
Following  summation  of  the  external  forces  and  solution  of  the  internal  (cable)  forces,  the 
solution  function  computes  the  generalised  accelerations  and  velocities  (it,  r)  at  the  current 
state.  This  point  solution  is  passed  to  the  integration  function  and  the  simulation  loop 
continues. 

It  is  also  possible  to  calculate  a  linear  model  by  numerical  approximation  of  the  Jaco- 
bians  yuM  and  from  the  nonlinear  system. 
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This  model  will  have  the  form 


u  =  [v«m]u  +  (29) 

and  can  be  used  for  an  alternative  linear  simulation  about  the  trim  state.  Another  use  is 
in  various  linear  system  analyses,  such  as  the  determination  of  the  natural  modes,  which 
will  be  discussed  in  the  following  section. 
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4  Simple  Pendulum  System  Dynamics 

The  first  phase  of  the  analysis  involved  an  investigation  of  the  dynamics  of  pendulum- 
type  systems.  For  this  purpose,  the  full  helicopter  slung-load  simulation  model  developed 
in  Section  3  was  constrained  so  as  to  approximate  a  two-body  pendulum  system.  Fur¬ 
thermore,  the  aerodynamic  effects  of  both  helicopter  and  load  were  excluded  from  the 
model. 

4.1  System  Properties 

The  simulation  was  validated  for  both  free  and  constrained  models  using  the  system 
given  in  Reference  [11]  for  a  CH-53D  helicopter  with  a  single  slung  load.  The  CH-53D  is 
a  twin-turbine,  main  and  tail  rotor  transport  helicopter  with  mass  and  inertia  properties 
as  outlined  below. 

Mass  (lb)  Moments  of  Inertia  (slug.ft2) 

TUh  I xx  Iyy  Izz  Ixz 

35000  36100  191500  179200  14800 

The  slung  load  chosen  was  a  standard  military  container,  known  as  a  MILVAN,  which 
is  a  common  helicopter  cargo  used  in  many  commercial  and  military  operations.  The 
dimensions  of  a  MILVAN  container  are  20  x  20  x  8  ft  and  the  mass  typically  varies  from 
4000  lb  (empty)  to  20000  lb  (full).  In  the  Reference  cited  above,  examples  with  masses 
outside  this  range  were  also  checked  to  demonstrate  very  low  density  and  very  heavy  loads. 
The  moments  of  inertia  for  the  container  load  were  approximated  with  linear  functions  of 
its  mass  by 

Ixx  =  0.33  *  mi  Iyy  =  Izz  =  1.20  *  mi  (30) 

where  the  moments  of  inertia  are  in  slug.ft2  and  the  mass  is  in  lb. 


4.2  Analysis  of  System  Modes 

In  order  to  validate  the  simulation  developed,  comparisons  were  made  against  those 
previously  reported  in  a  numerical  example  given  in  Reference  [11].  In  addition,  analyt¬ 
ical  results  were  calculated  for  the  modes  of  similar  pendulum  systems.  The  governing 
equations  for  these  systems  are  detailed  in  Appendix  B,  along  with  several  simplifying 
approximations  and  solutions  to  each  in  terms  of  their  natural  frequencies. 

For  the  example  cited,  the  mass  ratio  is  mi/mh  =  0.05.  Using  Equation  (30),  this 
yields  the  following  load  properties, 

Mass  (lb)  Moments  of  Inertia  (slug.ft2) 

mi  I xx  Iyy  Izz  Ixz 

1750  577.5  2100  2100  0 

The  sling  configuration  used  consisted  of  a  single  pendant  suspension  and  bridle,  as 
illustrated  in  Figure  4.  The  total  sling  length  L  between  helicopter  attachment  point  and 
the  load  eg  is  25  ft,  and  the  pendant-to-sling  length  ratio  i/L  is  0.6. 
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Figure  4:  Two-body  Pendulum  Simulation  System 


For  the  system  examined,  there  are  two  oscillatory  modes  in  both  longitudinal  and 
lateral  axes.  Essentially,  the  low  frequency  mode  is  associated  with  the  pendulous  motion 
of  the  load  along  the  total  sling  length.  The  higher  frequency  mode  is  associated  with  the 
coupled  pitching  (or  rolling)  motion  of  the  load  and  bridle  and  the  pendant  suspension. 
Values  for  the  natural  frequencies  obtained  from  the  analytical  formulation  and  the  sim¬ 
ulation  are  listed  in  Table  1.  The  first  column  of  analytical  results  was  calculated  from 
Equation  (B16).  The  second  and  third  columns  pertain  to  the  models  generated  by  the 
simulation  code  HSLSIM  developed  in  the  current  study.  Of  these,  one  was  constrained 
so  as  to  approximate  the  pendulum  system  used  in  the  analytical  derivation,  which  for 
the  longitudinal  case  included  constraints  in  translation  along  both  y  and  z  axes,  and 
constraint  in  rotation  about  the  y  axis.  The  other  model  was  free  to  move  in  all  axes,  pro¬ 
viding  a  closer  approximation  to  the  real  system.  The  fourth  column  lists  the  frequencies 
obtained  in  the  previous  analysis  using  the  simulation  code  EOMPROG. 


Table  1:  Pendulum  Mode  Frequencies  (rad/s)  for  CH-47B-MILVAN  System  in  Hover 


Mode 

Analytical 
(Eq.  (B16)) 

HSLSIM 

Constrained 

HSLSIM 

Free 

EOMPROG 
(Ref.  [11]) 

Longitudinal  (mode  I) 

1.12 

1.12 

1.15 

1.14 

Longitudinal  ( mode  II) 

3.86 

3.86 

3.86 

3.84 

Lateral-roll  (mode  I) 

1.15 

1.16 

1.25 

1.24 

Lateral-roll  (mode  II) 

7.17 

7.17 

7.18 

7.14 

Agreement  between  the  analytical  results  and  the  constrained  model  is  very  good. 
However,  there  are  some  small  discrepancies  between  those  first  two  sets  and  the  last 
two  sets  —  the  unconstrained  models.  This  can  be  explained  by  an  additional  coupling 
effect  in  the  unconstrained  models,  as  the  sling  force  at  the  attachment  point  produces 
a  moment  about  the  helicopter  eg.  The  effect  is  more  prominent  in  the  lateral  case, 
since  the  helicopter  moment  of  inertia  is  relatively  low  in  that  axis.  Differences  between 
the  unconstrained  model  generated  in  HSLSIM  and  EOMPROG  are  understandable,  as 
they  were  generated  by  two  quite  different  approaches.  The  simulation  code  HSLSIM 
incorporates  a  full  nonlinear  representation  of  the  helicopter  slung-load  system,  which  was 
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linearised  numerically  about  the  trim  state  to  obtain  a  Jacobian  matrix  for  model  analysis 
( see  the  following  section).  Using  this  Jacobian  matrix  approach  would  therefore  incur 
errors  from  the  numerical  approximation.  EOMPROG,  on  the  other  hand,  is  based  on 
an  explicitly  linear  small-perturbation  formulation  and,  consequently,  errors  would  arise 
from  such  simplification  as  the  small  angle  assumptions  and  the  exclusion  of  higher  order 
terms. 

It  is  also  interesting  to  see  how  the  analytical  solution  compares  with  the  unconstrained 
simulation  model  HSLSIM  for  a  range  of  configurations.  Figures  5  and  6  show  the  variation 
in  natural  frequency  for  the  longitudinal  pendulum  modes  with  load-to-helicopter  mass 
ratio  mi/mh  and  pendant-to-sling  length  ratio  l/L,  respectively.  Clearly,  there  is  quite 
a  large  difference  in  the  higher  frequency  mode  over  the  range  of  IjL  with  a  smaller 
deviation  in  the  lower  frequency  mode.  The  analytical  approximation  is  most  accurate 
at  a  length  ratio  of  0.6,  but  outside  this  differs  dramatically.  At  a  length  ratio  of  0.1, 
the  error  is  close  to  50%  of  the  frequency  predicted  by  HSLSIM.  At  a  length  ratio  of  0.9, 
the  error  is  approximately  120%.  There  is  less  difference  in  the  modes  over  the  range 
of  mi/mh.  For  the  configurations  examined,  the  largest  error  of  20%  occurred  in  the 
lower  frequency  mode  at  a  mass  ratio  of  1.0.  This  error  improved  as  the  mass  ratio  was 
decreased,  reaching  a  level  of  5%  at  a  ratio  of  0.1.  Both  of  these  discrepancies  are  due 
to  the  simplification  of  the  analytical  model,  which  prohibits  any  pitching  motion  about 
the  helicopter  eg.  Therefore,  as  a  result  of  this  test,  one  can  conclude  that  the  analytical 
approximation  given  in  Appendix  B  is  only  appropriate  for  low  load-to-helicopter  mass 
ratios  and  moderate  pendant-to-sling  length  ratios  (0.6  for  the  system  examined). 


4.3  Longitudinal  and  Lateral  Simulations 

Following  linear  analysis  of  the  modes,  several  simulations  were  run  so  as  to  gain  a 
better  picture  of  the  inherent  dynamic  behaviour  of  the  slung-load  system.  Simulations  of 
both  longitudinal  and  lateral  subsystems  were  conducted.  From  the  longitudinal  simula¬ 
tions,  two  cases  have  been  selected  for  inclusion  here  to  illustrate  the  effect  of  heavy  loads 
on  the  helicopter  response.  From  the  lateral  simulations,  two  cases  axe  included,  which 
demonstrate  the  instability  in  yaw  for  a  multiple-load  system. 

For  the  longitudinal  simulations,  the  helicopter-load  configuration  was  the  same  as  that 
used  in  the  linear  analysis,  i.e.  a  CH-53D  helicopter  and  MILVAN  load  attached  by  a  single 
pendant  suspension.  It  was  also  constrained  in  the  lateral-directional  axes,  as  before.  The 
slung  load  was  initially  offset  from  its  static  equilibrium  position  at  hover,  with  the  sling 
cable  set  at  -30°  in  pitch  and  the  load  at  -15°.  It  was  then  allowed  free  response  over  the 
10  s  duration.  The  time  history  response  plots  for  two  different  load  configurations,  with 
mass  ratios  of  0.05  and  1.0,  are  presented  in  Figures  7  and  8,  respectively.  In  the  first 
set  of  plots,  it  can  be  seen  that  the  load,  with  a  small  mass  compared  to  the  helicopter, 
has  very  little  effect  on  the  helicopter  motion,  as  expected.  Both  low  and  high  pendulum 
natural  frequencies  (0.18  and  0.61  Hz)  of  the  load  can  be  readily  identified  in  the  traces 
of  forward  velocity  and  pitch  rate,  respectively.  The  normal  acceleration  and  cable  force 
oscillate  at  approximately  twice  the  high  pendulum  natural  frequency  and  are  highly 
correlated,  with  peaks  at  the  pitch  rate  extremes  (0.36  Hz).  For  the  cable  force,  these 
peaks  have  magnitudes  which  are  roughly  1.5  times  the  static  load.  The  longitudinal 
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acceleration  is  also  highly  correlated  with  the  cable  angle.  In  the  second  set  of  plots,  it  is 
clear  that  the  load,  now  with  mass  equal  to  the  helicopter,  has  had  a  significant  effect  on 
the  helicopter  motion.  The  pendulum  natural  frequencies  are  also  less  well  defined  (0.30 
and  0.64  Hz).  Unlike  the  first  simulation,  in  this  case  the  helicopter  reaches  velocities  of 
similar  magnitude,  as  the  load  and  the  accelerations  in  forward  and  normal  directions  are 
out  of  phase  by  180°.  Once  again,  the  cable  force  oscillates  at  approximately  twice  the 
high  pendulum  natural  frequency,  with  peaks  near  the  pitch  rate  maxima  (0.33  Hz). 

For  the  lateral  simulations,  a  CH-47B  helicopter  and  two  box  containers  attached  by 
multiple  (bridle)  cables  were  used.  The  helicopter  was  constrained  in  the  longitudinal  and 
lateral  axes,  which  restricted  it  to  a  yawing/side  motion.  The  forward-attached  slung  load 
was  initially  offset  by  a  bank  angle  of  -30°  so  as  to  be  displaced  to  the  starboard  side  of 
the  helicopter,  and  the  aft-attached  slung  load  was  offset  by  an  angle  of  30°,  resulting  in 
a  displacement  to  the  port  side.  Once  again  the  system  was  allowed  free  response  over 
10  s.  The  time  history  response  plots  are  presented  in  Figures  9  and  10,  for  configurations 
with  total  mass  ratios  of  0.05  and  1.0-  From  these  plots,  the  symmetric  nature  of  one 
of  the  pendulum  modes  can  be  seen.  Specifically,  the  velocity  component  v,  bank  angle 
(j),  pitch  angle  6,  and  acceleration  component  v  for  each  load  is  ‘mirrored’  by  the  other. 
Slight  deviations  have  come  about  because  of  the  different  attachment  point  locations  with 
respect  to  the  helicopter  eg.  It  is  also  important  to  note  that  there  is  a  very  low  frequency 
response  in  yaw  for  both  helicopter  and  load,  which  manifests  itself  through  the  yaw  rate 
r  and  azimuth  angle  ip.  Although  these  two  traces  seem  to  be  diverging,  they  are  actually 
just  at  the  start  of  a  long  period  (<  0.01  Hz)  oscillation.  This  oscillation  has  an  amplitude 
of  approximately  160°  in  the  helicopter  azimuth  angle.  The  effect  of  the  loads  swinging  in 
this  manner  therefore  has  a  significant  influence  on  the  yaw  of  the  helicopter. 
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Figure  8: 


-  Helicopter  -  Slung-load 


CH-53D-MILVAN  System  Time  History  Response  —  Mass  Ratio  0.05 


CH-53D-MILVAN  System  Time  History  Response  —  Mass  Ratio  1.00 
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5  Helicopter  Slung-Load  Dynamics 

The  second  phase  of  the  analysis  examined  the  dynamics  of  the  full  helicopter  slung- 
load  simulation  model.  Both  the  CH-53D  and  the  CH-47B  helicopters  with  various  slung- 
load  configurations  were  considered.  These  models  were  free  in  all  axes  and  incorporated 
basic  helicopter  aerodynamics. 

The  CH-53D  helicopter  model  was  essentially  the  same  as  used  in  the  previous  section, 
but  with  the  inclusion  of  aerodynamic  effects.  Furthermore,  the  same  MILVAN  slung 
load  was  used  in  all  of  the  single-load  configurations  examined.  The  main  focus  of  the 
analysis,  however,  was  various  configurations  of  the  CH-47B  helicopter  with  slung  loads. 
The  CH-47B  is  a  twin-turbine,  tandem  rotor  transport  helicopter  with  multiple  sling 
attachment  points.  It  has  the  following  mass  and  inertia  properties: 

Mass  (lb)  Moments  of  Inertia  (slug.ft2) 

Trih  I xx  lyy  Izz  Xr  z 

33000  34000  202500  191000  14900 


5.1  Unloaded  Helicopter  Dynamic  Modes 

The  aerodynamic  model  used  for  both  helicopters  was  based  on  simple  linear  state- 
space  formulation  with  six  degrees  of  freedom  and  four  control  inputs.  For  the  CH-53D, 
the  stability  and  control  derivatives,  as  well  as  the  corresponding  trim  conditions,  were 
obtained  from  Reference  [35].  Several  small  modifications  were  also  made  as  applied  in 
Reference  [11]  in  order  to  perform  a  valid  comparison.  For  the  CH-47B,  the  derivatives 
and  trim  conditions  were  obtained  from  Reference  [27].  Since  the  aim  of  this  analysis 
was  simply  to  determine  the  modes  of  oscillation  for  the  helicopter  slung-load  system,  the 
Automatic  Flight  Control  System  (AFCS)  was  not  implemented. 

In  the  above  references,  the  stability  and  control  derivatives  for  both  helicopter  models 
were  tabulated  for  a  range  in  trim  airspeed.  For  the  CH-53D,  the  derivatives  were  given 
for  speeds  of  0,  10,  20,  40,  60,  80,  100,  120,  and  140  KTAS.  Furthermore,  to  illustrate  the 
change  in  behaviour  of  the  unloaded  helicopter,  the  eigenvalues  were  calculated  at  each  of 
these  airspeeds  and  drawn  on  the  complex  plane  in  the  form  of  a  root  locus  plot. 

Validation  of  the  helicopter  model  developed  in  the  current  work  was  performed  by 
calculating  the  same  eigenvalues  at  each  airspeed  and  then  comparing  against  those  re¬ 
ported  in  Reference  [11].  Figure  11  presents  the  eigenvalues  for  the  CH-53D  obtained  from 
both  analyses. 

Although  there  are  some  small  differences,  the  eigenvalues  are  generally  quite  close. 
As  in  the  previous  section,  those  differences  can  be  explained  by  the  different  linearisation 
approaches.  There  was  also  some  inherent  error  in  transcribing  the  published  results  to 
another  graph  by  hand.  The  behaviour  of  the  CH-53D  is  fully  analysed  in  Reference  [11] 
and  will  not  be  discussed  here.  Instead,  the  behaviour  of  the  CH-47B  over  a  range  of 
airspeeds  will  be  examined.  In  Reference  [27],  the  stability  and  control  derivatives  were 
given  for  speeds  of  0.1,  20,  40,  60,  80,  100,  120,  and  130  KTAS.  Again,  the  eigenvalues 
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Figure  11:  CH-53D  System  Eigenvalues  —  Variation  with  Trim  Airspeed  (in  KTAS) 


were  calculated  at  each  of  these  airspeeds  and  drawn  on  the  complex  plane,  as  shown  in 
Figure  12. 

The  gridlines  in  each  of  these  graphs  represent  lines  of  constant  frequency  and  constant 
damping  factor.  The  curved  lines  define  lines  of  constant  frequency,  which  increase  in  value 
from  zero  as  they  radiate  out  from  the  origin.  The  straight  lines  define  lines  of  constant 
damping  factor.  They  decrease  in  a  clockwise  direction  from  a  value  of  1.0  (along  the 
negative  real  axis)  to  a  value  of  -1.0  (along  the  posisive  real  axis).  The  positive  imaginary 
axis  represents  zero  damping.  Any  modes  on  the  positive  right-hand  side  of  the  real  axis 
thus  have  a  negative  damping  factor  and  are  therefore  unstable.  Eigenvalues  with  non-zero 
imaginary  components  represent  oscillatory  modes  with  a  certain  frequency  and  damping 
factor,  while  those  that  lie  along  the  real  axis  represent  pure  damping  modes. 

The  CH-47B,  like  most  conventional  helicopters,  has  essentially  three  longitudinal 
modes  and  three  lateral  modes.  However,  it  should  be  noted  that  these  modes  change 
markedly  with  airspeed  and  tend  to  be  highly  coupled.  The  modes  consist  of: 

•  A  low-frequency,  marginally  unstable  phugoid  mode  in  hover,  which  decreases  slightly 
in  damping  as  the  airspeed  increases  to  20  KTAS.  The  frequency  then  decreases  and 
becomes  stable  at  around  40  KTAS.  Between  40  and  60  KTAS,  the  frequency  in¬ 
creases  before  decreasing  again  through  to  130  KTAS.  There  is  a  distinct  kink  in 
the  root  locus  at  60  KTAS,  mainly  due  to  rapid  changes  in  the  influential  stability 
derivatives  in  this  airspeed  regime.  Above  80  KTAS,  the  mode  takes  on  a  signifi¬ 
cant  lateral  component,  or  side-motion.  At  130  KTAS,  the  mode  comprises  all  three 
translational  u,  v,  and  w  velocity  components  equally. 

•  A  higher  frequency,  marginally  unstable  lateral  phugoid  mode,  which  also  decreases 
slightly  in  damping  as  the  airspeed  increases  to  20  KTAS.  From  20  to  130  KTAS, 
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Figure  12:  CH-47B  System  Eigenvalues  —  Variation  with  Trim  Airspeed  (in  KTAS) 

the  frequency  increases  uniformly  and  the  damping  remains  relatively  constant.  In 
hover,  the  lateral  phugoid  mode  resembles  the  low-frequency  ‘swinging’  motion  of 
the  longitudinal  phugoid  mode.  As  the  airspeed  increases,  however,  the  mode  looks 
much  more  like  a  dutch-roll  mode,  typical  of  fixed-wing  aircraft. 

•  A  stable,  pure  damping  mode  in  pitch,  or  pitch  subsidence,  which  initially  decreases 
and  then  increases  in  damping  with  an  increase  in  airspeed.  This  mode  actually  has 
a  larger  roll-rate  component  at  hover,  then  as  the  airspeed  increases  the  pitch-rate 
component  becomes  more  dominant,  reaching  a  peak  at  40  KTAS  before  decreasing. 
Above  60  KTAS,  the  mode  has  a  large  w  velocity  component. 

•  Another  stable,  pure  damping  mode  in  roll,  or  roll  subsidence ,  which  initially  in¬ 
creases  and  then  decreases  in  damping  with  an  increase  in  airspeed.  Complimentary 
to  the  previous  mode,  this  mode  has  a  larger  pitch-rate  component  at  hover,  and  as 
the  airspeed  increases,  the  roll-rate  component  becomes  more  dominant.  There  is 
little  change  in  the  mode  shape  thereafter. 

•  A  stable,  pure  damping  heave  subsidence  mode,  which  generally  decreases  in  damp¬ 
ing  with  an  increase  in  airspeed.  The  u  and  v  velocity  components  increase  and  the  w 
velocity  component  decreases  from  hover  to  approximately  30  KTAS.  At  this  point, 
the  u  and  v  velocity  components  both  decrease  to  near-zero  values.  The  w  velocity 
decreases  by  only  a  small  amount  before  remaining  relatively  constant  through  to 
130  KTAS. 

•  An  initially  stable  pure  damping  mode  in  yaw,  which  decreases  in  damping  with  an 
increase  in  airspeed,  becoming  statically  unstable  at  approximately  32  KTAS.  The 
damping  continues  to  decrease  rapidly  until  an  airspeed  of  40  KTAS,  when  it  settles 
to  a  constant  value.  The  v  velocity  component  varies  throughout  the  speed  range, 
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while  the  u  any  w  components  both  increase  greatly  at  the  crossover  airspeed.  The 
pitch  rate  component  also  increases  at  32  KTAS,  becoming  the  dominant  part  of 
the  motion.  The  yaw  rate  component,  however,  actually  decreases  from  its  value  at 
hover  and  remains  small  by  comparison. 


The  eigenvalues  representing  each  of  these  modes  at  airspeeds  of  0.1  KTAS  (hover) 
and  130  KTAS  are  tabulated  below.  The  corresponding  natural  frequency  and  damping 
factor  are  also  listed  for  each  mode. 


Table  2:  Natural  Frequency  and  Damping  for  CH-4 7B  at  0.1  KTAS 


Mode 

Eigenvalue 

Natural 

Frequency  (rad/s) 

Damping 

Factor 

Longitudinal  phugoid 

0.1099  ±  0.5026* 

0.5145 

-0.2137 

Pitch  subsidence 

-1.4853 

1.4853 

1.0000 

Heave  subsidence 

-0.3003 

0.3003 

1.0000 

Lateral  phugoid 

0.0453  ±  0.4829i 

0.4850 

-0.0934 

Roll  subsidence 

-1.3396 

1.3396 

1.0000 

Yaw  mode 

-0.0766 

0.0766 

1.0000 

Table  3:  Natural  Frequency  and  Damping  for  CH-47B  at  130.0  KTAS 


Mode 

Eigenvalue 

Natural 

Frequency  (rad/s) 

Damping 

Factor 

Longitudinal  phugoid 

-0.0619  ±  0.1534* 

0.1654 

0.3744 

Pitch  subsidence 

-2.9048 

2.9048 

1.0000 

Heave  subsidence 

-0.0144 

0.0144 

1.0000 

Lateral  phugoid 

0.0610  ±  0.8754* 

0.8775 

-0.0695 

Roll  subsidence 

-1.2224 

1.2224 

1.0000 

Yaw  mode 

0.6008 

0.6008 

1.0000 

5.2  Analysis  of  Helicopter  Slung-Load  System  Modes 

As  with  the  unloaded  helicopter  model,  the  results  obtained  for  the  combined  helicopter 
and  load  were  compared  with  those  documented  in  Reference  [11].  The  eigenvalues  for  a 
range  of  configurations  of  the  CH-53D  in  hover  with  a  MILVAN  load  were  examined.  In 
particular,  the  load-to-helicopter  mass  ratio  mi/m^  was  varied  between  0.1  and  0.6.  The 
load  was  slung  using  multiple  cables  (sling  legs)  and  a  single  helicopter  attachment  point, 
and  the  distance  between  the  attachment  point  and  load  eg  was  25  ft.  Figures  13  and  14 
illustrate  the  root  loci  for  this  system  in  both  longitudinal  and  lateral  axes. 

The  additional  eigenvalues  in  these  diagrams  represent  the  modes  associated  with  the 
pendulous  motion  of  the  load  and  will  be  discussed  later  in  this  section.  Note  also  that 
the  heave  subsidence  and  yawing  modes  were  not  included  in  the  previous  results  and 
are  therefore  not  shown  here.  Nonetheless,  all  of  the  remaining  modes  agree  quite  well, 
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Figure  13:  CH-53D  with  Single  Slung  Load  Longitudinal  Eigenvalues  at  Hover  —  Variation 
with  Load-to-Helicopter  Mass  Ratio 


Figure  14:  CH-53D  with  Single  Slung  Load  Lateral  Eigenvalues  at  Hover  —  Variation  with 
Load-to-Helicopter  Mass  Ratio 
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particularly  in  the  lateral  axes.  Small  discrepancies  have  come  about  due  to  the  same 
reasons  discussed  above. 

A  similar  numerical  analysis  was  conducted  with  the  CH-47B  in  hover  and  a  single 
MILVAN  load  attached  with  the  same  sling  configuration  as  the  CH-53D,  shown  in  Fig¬ 
ure  15. 


Figure  15:  Multiple  Cable  Sling  Configuration 

The  load-to-helicopter  mass  ratio  in  this  case  was  varied  from  0.1  to  1.0.  Figures  16 
and  17  illustrate  the  change  in  system  behaviour  with  mass  ratio. 

The  helicopter  slung-load  system  comprises  the  same  typical  helicopter  modes  intro¬ 
duced  previously,  as  well  as  an  additional  pendulum  mode  in  each  axis.  These  pendulum 
modes  can  be  attributed  to  the  simple  swinging  motion  of  the  load  beneath  the  helicopter. 
For  low  load  weights,  the  modes  are  largely  decoupled  from  the  helicopter  motion  and 
approach  that  of  a  similar  pendulum  system,  with  frequencies  as  given  in  Appendix  B, 
Equation  (B23).  As  the  load  weight  is  increased,  it  has  more  influence  on  the  motion  of 
the  helicopter.  The  frequencies  approach  that  of  a  double-pendulum  system,  as  in  Equa¬ 
tion  (B21),  in  the  appendix2.  Furthermore,  the  damping  for  both  lateral  and  longitudinal 
pendulum  modes  generally  increases,  producing  a  stabilising  effect. 

In  contrast,  the  effect  of  increasing  the  load  mass  on  the  other  modes  is  destabilising. 
The  damping  of  the  pitch,  heave,  and  roll  modes  generally  decreases;  the  yaw  mode  changes 
very  little.  The  effect  on  the  unstable  phugoid  modes  is  perhaps  more  critical.  From  the 
diagrams,  the  longitudinal  phugoid  mode  mainly  changes  in  frequency  —  only  decreasing 
by  a  small  amount.  The  lateral  phugoid  mode,  however,  undergoes  a  significant  decrease 
in  damping  as  the  load  weight  increases.  Since  this  mode  is  already  unstable,  this  would 
result  in  a  more  rapid  divergence  from  trim.  In  summary,  the  net  effect  of  increasing  the 
load  weight  would  therefore  be  a  marginal  stabilisation  of  the  pendulum  modes  and  a  more 
significant  destabilisation  of  the  phugoid  modes. 

Another  numerical  experiment  with  the  CH-47B  in  hover  and  a  single  MILVAN  load 
attached  with  a  range  of  sling  configurations  was  also  conducted.  In  this  test,  the  load  and 
bridle  arrangement  was  attached  by  a  single  pendant  cable  to  the  helicopter,  as  shown  in 
Figure  18. 

The  ratio  of  the  pendant  length  to  total  sling  length,  i/L ,  was  varied  from  0.1  to  0.9. 

2  A  more  accurate  estimate  which  incorporates  the  helicopter’s  moment  of  inertia  can  be  obtained  by 
using  Equation  (B15). 
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Figure  16:  CH-47B  with  Single  Slung  Load  Longitudinal  Eigenvalues  at  0.1  KTAS  — 
Variation  with  Load-to-Helicopter  Mass  Ratio 


Re(s) 

Figure  17:  CH-47B  with  Single  Slung  Load  Lateral  Eigenvalues  at  0.1  KTAS  —  Variation 
with  Load-to-Helicopter  Mass  Ratio 
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Figure  18:  Single  Cable  Sling  Configuration 


The  total  sling  length  (the  distance  between  the  helicopter  attachment  point  and  load  eg) 
remained  constant  at  25  ft  and  the  mass  ratio  was  0.05.  Root  locii  for  the  longitudinal 
and  lateral  modes  are  illustrated  in  Figures  19  and  20,  respectively. 

In  this  configuration,  there  are  two  pendulum  modes  in  each  axis,  which  reflect  the 
extra  two  degrees  of  freedom  in  load  pitch  and  load  roll.  The  first  low-frequency  mode 
can  be  associated  with  the  gross  swinging  motion  of  the  load-bridle  subsystem  beneath 
the  helicopter.  The  second  high-frequency  mode  can  be  associated  with  the  pitching  or 
rolling  motion  of  the  load  about  its  eg.  Since  the  mass  ratio  in  this  case  is  quite  low,  the 
natural  frequencies  of  oscillation  can  be  estimated  quite  accurately  using  Equation  (B15) 
in  Appendix  B. 

As  the  length  ratio  increases,  this  second  pendulum  mode  increases  in  frequency  quite 
dramatically.  For  this  particular  load,  the  frequency  of  the  lateral  mode  is  approximately 
double  that  of  the  longitudinal  mode  due  to  the  smaller  moment  of  inertia  about  the  roll 
axis.  All  other  modes,  including  the  low  frequency  pendulum  modes,  do  not  change  by 
any  significant  amount  with  an  increase  in  the  length  ratio. 

The  last  component  in  the  analysis  of  the  system  modes  was  concerned  with  multiple¬ 
load  configurations.  To  this  end,  the  dynamic  characteristics  of  the  CH-47B  with  multiple, 
individually  slung  loads  was  examined.  It  could  be  reasonably  assumed  that  the  addition 
of  extra  loads  would  simply  add  more  pendulum  modes  to  the  system.  However,  it  is 
interesting  to  see  the  placement  of  those  modes,  as  well  as  any  effect  on  the  other  modes. 
Figure  21  illustrates  the  variation  in  the  system  modes  for  zero,  one,  two,  and  three 
individually  attached  loads.  The  parameter  n  denotes  the  number  of  bodies,  including  the 
helicopter,  in  the  system. 

For  this  analysis,  simple  box-type  loads  were  used.  As  with  the  MILVAN,  the  moments 
of  inertia  for  each  load  were  approximated  with  linear  functions  of  their  mass  by 

I XX  =  Iyy  =  Izz  =  0.33  *  mi  (31) 

The  weight  of  each  load  was  the  same,  and  for  each  configuration,  the  total  load-to- 
helicopter  mass  ratio  was  kept  constant  at  1.0.  The  loads  were  each  slung  using  multiple 
cables  (sling  legs)  attached  to  individual  helicopter  hook  points,  and  the  distance  between 
each  attachment  point  and  the  corresponding  load  eg  was  25  ft. 
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Figure  19:  CH-47B  with  Single  Slung  Load  Longitudinal  Eigenvalues  at  0.1  KTAS  — 
Variation  with  P endant-to- Sling  Length  Ratio 
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Figure  20:  CH-47B  with  Single  Slung  Load  Lateral  Eigenvalues  at  0.1  KTAS  —  Variation 
with  Pendant-to-Sling  Length  Ratio 
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Figure  21:  CH-47B  with  Multiple  Slung  Loads  Eigenvalues  at  0.1  KTAS  —  Variation  with 
Number  of  Loads 

For  the  unloaded  helicopter,  the  modes  consist  of  two  oscillatory  modes  (longitudi¬ 
nal  and  lateral  phugoid)  and  four  damping  modes  (including  pitch,  roll,  heave,  and  yaw 
subsidence).  Recall  from  the  previous  section  that  the  pitch  and  roll  modes  are  actually 
swapped  around  in  hover.  The  single-load  (centre  hook)  configuration  has  two  additional 
modes  (one  longitudinal  and  one  lateral  pendulum  mode),  both  of  which  are  stable.  How¬ 
ever,  the  addition  of  one  or  more  loads,  or  a  number  of  loads  further  destabilises  the  lateral 
phugoid  mode.  It  also  significantly  decreases  the  damping  in  pitch,  roll,  and  heave  modes 
and,  to  a  much  lesser  extent,  the  yaw  mode.  The  dual-load  (forward  and  aft  hooks)  config¬ 
uration  has  those,  plus  two  more  oscillatory  modes,  which  are  all  highly  coupled.  The  first 
two  pendulum  modes,  which  have  approximately  the  same  locations  as  in  the  single-load 
case,  describe  longitudinal  and  lateral  motions  in  which  the  loads  swing  in  phase.  The 
next  two  pendulum  modes  are  very  lightly  damped  and  describe  motions  in  which  the 
loads  swing  out  of  phase.  They  have  little  effect  on  the  helicopter  motion.  The  triple-load 
(all  hooks)  configuration  also  has  the  typical  single-load  style  modes,  plus  four  additional 
pendulum  modes.  They  describe  two  longitudinal  and  two  lateral  motions  in  which  the 
loads  swing  in  various  combinations. 
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5.3  Flight  Simulation 

Further  to  the  modal  analysis,  a  number  of  simulations  were  run  in  order  to  demon¬ 
strate  the  typical  behaviour  of  helicopter  slung-load  systems.  The  first  of  these  demon¬ 
strates  the  longitudinal  response  of  a  simple  single-load  configuration.  In  this  case,  a 
MILVAN  load  was  slung  with  a  single  pendant  cable  to  a  CH-47B  helicopter,  as  in  the 
configuration  shown  previously  (Figure  18).  The  ratio  of  the  pendant  length  to  total  sling 
length  was  set  to  0.6,  and  the  total  sling  length  to  25  ft.  The  load-to-helicopter  mass  ratio 
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was  0.05. 

Figure  22  shows  time-histories  of  the  body-axis  velocities,  pitch  rate  and  angle,  cable 
angle  and  tension,  and  control  inputs.  Figure  23  illustrates  the  motion  of  the  helicopter 
and  slung  load  through  the  duration  of  the  simulation. 

In  the  time-history  plots,  u  and  w  are  the  velocities  in  x  and  2  directions  of  the  body 
axes  respectively,  q  is  the  pitch  rate,  0  is  the  pitch  angle,  and  6C  is  the  cable  angle  displaced 
from  the  vertical  position.  The  cable  tension  force,  nondimensionalised  by  the  load  weight, 
is  denoted  by  fc.  The  longitudinal  cyclic  and  collective  control  inputs  are  represented  by 
5b  and  5C,  respectively. 

For  this  10  second  simulation,  the  helicopter  and  slung  load  were  given  an  initial 
forward  velocity  of  80  ft/s  and  the  load  was  displaced  from  its  static  equilibrium  position 
by  30°.  The  control  inputs  were  generated  so  as  to  maintain  a  nose- level  attitude  in  the 
helicopter.  Both  primary  and  secondary  pendulum  modes  can  be  readily  identified  in  the 
response  of  the  load,  most  notably  in  the  cable  angle.  In  general,  the  system  behaved  as 
expected,  with  the  load  having  little  effect  on  the  behaviour  of  the  helicopter.  It  is  also 
worth  noting'that  the  peaks  in  the  cable  tension  correspond  to  the  extremes  of  the  pitch 
rates  and  have  a  maximum  value  of  up  to  1.5  times  their  static  load. 

The  second  simulation  demonstrates  the  combined  longitudinal  and  lateral  response  of 
the  same  single-load  configuration.  However,  for  this  case,  the  mass  ratio  was  increased 
to  1.0.  The  duration  was  20  seconds,  although  only  the  first  10  seconds  of  the  simulation 
are  shown  in  Figure  24.  Variables  shown  are  the  body-axis  velocities  and  accelerations, 
the  angular  displacements  and  rates,  and  the  cable  angles  and  tension.  The  flight  path  of 
the  helicopter  and  slung  load  are  depicted  in  Figure  25. 

In  these  plots,  v  is  the  velocity  in  the  y  direction  of  the  body  axes,  p  and  r  are  the  roll 
and  yaw  rates,  is  the  bank  angle,  and  <pc  is  the  cable  angle  displaced  from  the  vertical 
position.  The  variables  u,  v,  and  u  represent  the  accelerations  in  body  axes. 

Unlike  the  previous  simulation,  the  initial  state  was  hover.  The  load  was  displaced 
from  its  static  equilibrium  position  by  30°  in  pitch  and  15°  in  bank.  In  order  to  observe 
the  influence  of  the  load  on  the  helicopter,  the  controls  were  held  fixed  (at  zero)  for  the 
entire  20  second  simulation.  Both  primary  and  secondary  longitudinal  pendulum  modes 
are  visible  in  the  cable  pitch  angle,  and  both  lateral  modes  are  also  visible  in  the  cable  bank 
angle.  These  also  induce  a  yawing  motion,  seen  in  the  yaw  rate,  due  to  lateral-directional 
coupling,  and  although  it  may  seem  so,  this  motion  is  not  oscillatory.  The  rolling  modes 
are  generally  of  higher  frequency  than  the  pitching  modes,  given  their  lower  moment  of 
inertia.  However,  due  to  the  more  unstable  lateral  phugoid  mode  of  the  helicopter,  the 
rolling  motion  of  the  load  has  a  greater  consequence  on  the  resultant  motion  of  the  system. 
This  is  what  causes  the  rapid  divergence  in  the  velocity  v.  Again,  the  cable  tension  reached 
a  peak  of  nearly  1.5  times  the  static  load. 

The  last  simulation  that  has  been  included  in  this  report  demonstrates  the  full  response 
of  a  multiple-load  configuration.  Three  box-shaped  containers  with  an  equal  mass  ratio 
of  0.33  were  attached  to  the  three  hook  points  under  the  CH-47B.  They  were  each  slung 
using  four  cables  (Figure  15),  and  the  length  from  the  helicopter  attachment  points  to 
each  load  eg  was  set  to  15,  20,  and  25  ft. 
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Figure  22:  Longitudinal  Dynamic  Response  Simulation  of  CH-/f  7B  with  Single  Slung  Load 
—  Time  Histories  (mi/mh  —  0.05,) 


Figure  23:  Longitudinal  Dynamic  Response  Simulation  of  CH-4  7 B  with  Single  Slung  Load 
—  Simulation  Frames  ( mi/m h  =  0.05,) 
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Figure  25:  Coupled  Dynamic  Response  Simulation  of  CH-f7B  with  Single  Slung  Load  — 
Simulation  Frames  (mi/mh  =  l.Oj 


Figure  26:  Coupled  Dynamic  Response  Simulation  of  CH-f7B  with  Three  Slung  Loads  — 
Simulation  Frames  (mi/mh  =  0.33,) 
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Figures  26  and  27  contain  the  simulation  time  histories  and  helicopter  slung-load  flight 
path  respectively.  Again,  only  the  first  10  seconds  of  the  time  histories  are  shown. 

As  before,  the  initial  state  for  the  helicopter  was  hover.  All  of  the  loads  were  displaced 
from  their  equilibrium  positions  in  both  pitch  and  bank  by  arbitrary  angles  from  0°  to  30°. 
The  controls  were  also  held  fixed  for  the  duration  of  the  simulation,  which  was  just  over 
20  seconds.  From  the  first  set  of  plots,  the  fundamental  frequency  for  each  of  the  slung 
loads  is  visible.  The  slight  difference  in  their  frequency,  due  to  the  different  sling  lengths, 
is  perhaps  best  seen  in  the  pitch  rate  and  angle.  There  is  also  a  significant  coupling  effect, 
which  induces  a  yawing  motion.  Unstable  longitudinal  and  lateral  modes  again  cause 
divergence  in  the  velocities  u  and  v. 
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6  Concluding  Remarks 

A  simulation  model  was  produced  using  the  equations  of  motion  for  general  slung-load 
systems  developed  by  Cicolani  and  Kanning  [29].  The  formulation  used  is  based  on  the 
Newton- Euler  equations  written  in  terms  of  generalised  coordinates  and  can  be  readily 
adapted  to  systems  with  either  elastic  or  inelastic  suspension.  All  code  development, 
including  the  simulation  routine,  linear  analysis,  and  graphical  replay  tools,  was  done 
in  MATLAB.  Further  to  the  model  described,  the  simulation  has  also  been  extended  to 
incorporate  multiple-load  systems. 

Following  comparison  of  the  behaviour  of  a  simple  pendulum-type  system  against  both 
analytical  solutions  and  previously  published  results,  the  open-loop  characteristics  of  a 
full  helicopter  slung-load  system  was  examined.  The  results  presented  include  the  modal 
analysis  of  a  single-load  system  and  the  simulation  of  a  multiple-load  system.  It  was  found 
that  the  frequency  of  the  longitudinal  and  lateral  pendulum  modes  generally  increases  with 
the  load-to-helicopter  mass  ratio.  More  importantly,  the  lateral  phugoid  mode  becomes 
significantly  unstable  as  the  mass  ratio  is  increased.  For  the  simulation  demonstrated,  the 
dominant  natural  frequencies  were  identified. 
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Appendix  A:  Governing  System  Equations 


The  general  helicopter  slung- load  system  under  consideration  is  illustrated  in  Figure  Al. 
As  explained  in  Section  3,  the  system  consists  of  a  single  helicopter  supporting  one  or 
more  loads  by  means  of  some  suspension.  The  model  is  comprised  of  n  rigid  bodies, 
Bl,B2,...,Bn ,  with  m  links  supporting  a  single  force.  If  the  loads  all  have  either  multiple 
or  single  cable  suspensions,  then  m  =  n  -  1.  Furthermore,  if  the  suspensions  are  all  of 
single-cable  type,  as  in  the  example  shown,  there  are  m  cables,  C2,...,Cn. 


General  Transformation  Matrices 

For  scalar  representation  of  the  vector  cross  products,  the  general  form  of  the  skew- 
symmetric  matrix  S(Va)  is  defined  as 

5(14)  =  [  (V  x  ia)fl  (V  x  ja)Q  (V  x  ka)a  ] 

V az  Vay 
0  Vqx 

V ax  0 

where  the  columns  of  5(14)  represent  cross  products  of  V  with  the  axes  of  the  frame  JFa, 
each  referred  to  Ba. 

Using  this  formulation,  the  cross  products  representing  Coriolis  velocities  and  acceler¬ 
ations  can  be  written  as 


0 

Vaz 

~vay 
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and  the  centrifugal  accelerations  as 

(cj  x  <jj  x  R)q  =  S2(ua )  Ra  =  -S(uja)  S(Ra)  cja 


Defining  the  Euler-angle  transformation 


Tb,N 


cos  tpb  cos  6h  sin  ifib  cos  9b  -  sin  9b 

sin  (fib  cos  ifib  sin  9b  -  cos  (fib  sin  ifib  sin  (fib  sin  ifib  sin  9b  +  cos  (fib  cos  ifib  sin  (fib  cos  9b 
cos  (fib  cos  ifib  sin  9b  +  sin  <fib  sin  ifib  cos  (fib  sin  ifib  sin  9b  —  sin  (fib  cos  ifib  cos  (fib  cos  9b 


and  the  inverse 


TN,b  =  T£N 


the  angular  velocities  can  be  expressed  in  terms  of  the  Euler-angle  rates  via  the  transfor¬ 
mation 

wbb  =  Wbbab 


where 


with  inverse 


Wbb  = 


1  0  —  sin  9  b 

0  cos  (fib  sin  (fib  cos  9b 
0  —  sin  (fib  cos  (fib  cos  9b 


Wb;1  = 


1  sin  (fib  tan  9b  cos  (fib  tan  9b 
0  cos  (fib  —  sin  (fib 

0  sin  (fib  /  cos  9b  cos  (fib /  cos  9b 


Cable- Axes  (Single-Cable  Suspension) 

Cable-direction  vector: 

kcjx  =  ~~  (Rj*N  +  Tw,jRj*jj  -  -Rljv  -  TN,\Rl*aji) 
ij  o 

Cable  angles: 

0 
0 
1 

cos  (ficj  sin  9 cj 
—  sin  (ficj 

COS  (ficj  COS  9 cj 

The  cable  angles  are  obtained  from  the  cable  direction  vector 

(ficj  =  sin”1(-/ccjiv(2)) 

9Cj  =  sin-1  (  kcjN(l)  /  cos(fiCJ) 


kcjN  =  TNiCjkcjcj  = 


cj 


cos  9, 
0 

—  sin  9, 


cj 


sin  ficj  sin  9 cj 

COS  (ficj 

sin  (ficj  cos  9cj 


cos 


t>cj  sin  9cj 


sin<fi 


CJ 


COS  (ficj  COS  9cj 
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Velocity  of  each  cable: 


where 


V ajjcj  =  lj  kcjcj  +  lj  kcjcj 

lj  cos  4>cj  0  0 

=  o  -lj  0 

0  0  1 


kcjcj  =  ucjCj  x  kcjCj 


4>cj 

COS  C? cj 
@cj  Sia  4>cj 


Qcj 

4>cj 

.  ij  J 


'  0  ‘ 

X 

0 

1 

0cj  cos  (f>cj 

~4>cj 

o 


so  that  the  cable  angle  rates  can  be  derived  as 


4>cj  —  V ajjcj (2)  /  lj 

9cj  —  V  ajjcj  (1)  /  lj  cos  <pcj 


Simulation  Equations  for  an  Arbitrary  Number  of  Loads 

Cable  lengths: 

lj(r)  =  \R1*n  +  TNtiRl*aj\  -  Rj*N  -  TN,jRj*jj\ 


The  eg  position  of  each  load: 

Rj*i\'  —  R1*n  +  Tn,\RV  aj\  +  Rajjiy  —  T^^Rj  jj 

and  their  velocity: 

Vj*N  =  VrN  +  fNtlRl*ajx  +  VajjN  -  fN,jRj*jj 

=  V1*N  -  TNAS{Rl*aji)wh  +  VajjN  +  TN>jS{Rj*jj)wjj 

where 

VajjN  =  TN,cjV  ajjcj 


Hence, 


V  j*N  =  V1*N  +  TncjV  ajjcj  +  Ajtn+\W\\  +  Ajtn+jWjj 


and  inversely 


V ajjcj  =  —TCJjnV  1*v  +  Tcj,NV  j*N  +  Bjtn+iwli  +  BJtn+Jwjj 
with  the  elements 


Aj<n+i  =  —TN,iS(Rl*aji) 
Ajjn+j  =  TN,jS(Rj  jj) 
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Bj,n+ 1 
Bj,n+j 


Tcj^N  Ajtn+l 
~TcjitfAj>n+j 


Configuration  kinematics: 


v  =  Au 


[  VI *N  1 

V2*n 

Vn*N 

wl\ 

lo22 

wnn 

L 

0 

Tn,c2 

0 

0 


0  I 
0  I 


Tn,cti  | 
+ 


0  0  •••  0  ' 

r  vi*  i 

■^2,71+1  -'4-2,71+2  0 

Va22C2 

An,  71.+  I  0  AU:2n 

Vanricn 

I 

toll 

io22 

„ 

wnn 

u  =  A  lv 


'  vi*N  • 

'  I 

0 

0 

1  0 

0 

0  ■ 

r  vi*n  1 

V  a22C2 

~Tc2,N 

Tc2,N 

0 

1  -62,71+1 

-62, 71+2 

0 

V2*n 

V  anncn 

-Tcn,N 

0 

Tan,N 

1  -677,71+1 

0 

Bn  fin 

Vn*N 

— 

= 

— 

— 

-  - 

+ 

— 

—  — 

— 

ioli 

10I1 

io22 

0 

I 

io22 

wnn 

_ 

_ 

wnn 

Acceleration  vector: 

v  —  Au  +  Au 


Applied  forces  and  inertia  coupling: 

fo=fg  +  fa  —  X  —  DAu 
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FOIn 

'  migN  ' 

'  FAIn  " 

0 

0 

F02n 

m29  N 

FA2N 

0 

m2A2u 

FOtin 

mn9N 

FAritf 

0 

rrmAnu 

— 

= 

— 

+ 

— 

— 

— 

— 

— 

MOli 

0 

MAh 

5(wli)Jlwli 

0 

M022 

0 

MA22 

S(u}22)J2u22 

0 

M0nn 

0 

MAnn 

_  S(umn)Jnumn 

0 

L  - '<■  -1 

where 

U 

■  0 

”  _J 

0  ••• 

- -J  L-  \ 

0  1  0 

0 

0  ' 

0 

^2,2  •  •  • 

0  I  ^2,71+1 

1  • 

42,71+2  •  •  • 

0 

A  = 

0 

0 

1  .  : 

An,n  |  -'4.n,n+ 1 

-  + 

1 

0 

-4-7i,2n 

0 

0 

where  the  elements 

Aj,j 

=  TffyCj  S{wcjcj ) 

Aj>n+ 1  =  -TNtiS(wli)S{Rl*aji) 
Aj,n+j  =  TN,jS(wjj )S(Rj  jj) 


Suspension  forces: 

fc  —  Hr 


‘  FCIn  ' 
FC2n 

kc2 h  •  ■ • 

—kc2tf 

kcriN 

0 

FCtin 

0 

—kcriN 

MCh 

MC22 

C12i 

-m2 

£lni 

0 

_  MCnn  _ 

0 

-£nnn 

where 

£lji  =  ( Rl*aj  x  kcj) i  =  S{Rl*aj)iTitN  kcju 
Ujj  =  ( Rj*j  x  kcj)j  =  S{Rj*j)jTjyNkcjN 

yielding 
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DSTO-TR-1257 


Total  specific  force: 


sf  =  D  \fo  +  fc) 


'  SFlN  ' 
SF2n 

'  (F01N+FClN)/ml  ' 
(F02^r  +  FC2^)jm2 

SFritf 

( FOtin  +  FCriH)/mn 

SMh 

SM22 

Jl-x(M01i  +MCli) 
J2-1(M022  +MC22) 

_  SMnn 

Jn_1(M0nn  +  MCnn) 

Configuration  acceleration: 

u  =  A~1sf 


Configuration  position  vector: 


r  — 


r  ^  1 

\  1 

i?2^ 

V2*n 

Rn*N 

-  f 

Vn*N 

al 

J 

Wl^wlx 

a2 

W2^w22 

an 

Wn~lvmn 
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Appendix  B:  Analytical  Results  for 
Pendulum-Type  Systems 


There  are  a  number  of  publications,  including  Reference  [36],  which  list  the  natural  fre¬ 
quencies  and  mode  shapes  for  various  pendulum-type  systems.  However,  for  the  purpose 
of  this  work,  the  equations  of  motion  were  derived  and  the  resulting  natural  frequencies 
found  for  a  set  of  systems  that  can  be  used  to  approximate  a  helicopter  and  slung  load. 
The  generalised  system  is  analysed  first  and  then  several  approximations  for  both  single 
and  multiple  cable  configurations  are  examined. 


Generalised  Two-Body  system 

The  generalized  two-body  system  examined  in  this  appendix  consists  of  a  sliding  mass 
which  is  pinned  at  its  centre  of  gravity  and  a  swinging  mass  attached  below  it,  as  illustrated 
in  Figure  Bl.  The  supporting  body  of  mass  ms  and  inertia  Js  is  free  to  pitch  and  transverse 
along  the  x-axis,  while  the  suspended  body  or  load  of  mass  mi  and  inertia  J/  is  constrained 
just  to  swing  about  its  attachment  point. 


Figure  Bl:  Generalized  Two-body  Pendulum-type  System 

In  order  to  solve  for  the  natural  frequencies  of  the  system,  the  equations  of  motion 
must  first  be  stipulated.  For  simplicity,  any  perturbations  about  the  equilibrium  position 
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are  assumed  to  be  small.  As  a  result,  the  trigonometric  components  of  a  small  angle  9  can 
be  reduced  to  sin#  ~  9  and  cos 9  ~  1. 

The  sum  of  moments  about  the  cgs  of  the  support  and  load  are  respectively 

Jse  =  Qr  (Bl) 

Jl§i  —  -Pa(0i  -  9)  +  Qa 

and  the  balance  of  horizontal  forces  for  each  body  are 

msx  =  P6  +  Q  (B2) 

mixi  =  -P6  -  Q 


The  normal  component  P  of  the  interaction  force  can  be  approximated  by  the  mass- 
force  of  the  load,  whereas  the  tangential  component  Q  is  given  by  the  sum  of  moments 
about  the  supporting  body  (Eq.  (Bl)),  i.e. 


P  ss  mig 

Q  =  —0 

r 


(B3) 


Substituting  these  into  Equation  (Bl)  produces  the  first  equation  of  motion,  which  can 
be  rewritten  as 

Ji&i  =  -miga(9i  -  9)  +  —9 


+  mma9  —  J$i  —  mma9i  =  0 
r 

The  second  equation  of  motion  is  derived  by  subtracting  Equations  (B2)  to  give 
xi-x  =  ——(P9  +  Q)~  — ( P9  +  Q) 

mi  ms 

and  applying  the  geometric  relationship 

xi  —  x  =  r9  +  a9i 


(B4) 

(B5) 

(B6) 

(B7) 


This  produces  the  second  equation  of  motion,  written  as 

r9  +  a9i  =  — — ( mig9  +  —9) - —  ( mig9  +  —9) 

mi  r  ms  r 

Js 


where 


(r  4 - km)9  +  gkm9  +  a9i  =  0 

m/r 


km  =  (1  +  — ) 

m  q 


In  matrix  format,  Equations  (B5)  and  (B9)  can  be  expressed  as 


Jsd 

r 

M. b 

mirKr) 


-Jl  ■ 

'  9  ' 

1 

miga 

—miga 

‘  9  ' 

'  0  ‘ 

a 

A  . 

1 

gkm 

0 

.  Ql  . 

0  _ 

(B8) 

(B9) 

(BIO) 

(Bll) 
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Assuming  that  the  free  vibrations  are  harmonic  in  time  with  frequency  u,  the  equations 
of  motion  will  have  a  non-trivial  solution  if  the  following  determinant  is  equal  to  zero: 


_  +  rriiga  J/w2  -  rrnga 

~(r  +  +  gkm 


(B12) 


This  results  in  the  characteristic  equation 

Acj4  +  Bu2  +  C  =  0 


where 


A 

B 

C 


JcCL 


J, 


+  Ji(r  -\ - km ) 

r  mir 

-miga2  -  Jigkm  -  miga(r  H - —km) 

mir 

Tnig2akm 


(B13) 


(B14) 


Finding  the  roots  of  Equation  (B13)  provides  a  solution  for  the  natural  frequencies, 

2  =  -B  ± | g  -  4 AC\]»  (B15) 

2  A 

This  general  system  can  be  used  to  approximate  a  helicopter  slung-load  configuration, 
where  the  helicopter  is  free  to  pitch  about  its  eg  and  the  load  is  attached  using  a  multi-cable 
sling.  It  is  also  possible  to  further  simplify  the  system  to  represent  other  configurations 
as  demonstrated  in  the  following  sections. 


Single-cable  System 

In  order  to  simulate  a  single-cable  slung-load  configuration,  shown  in  Figure  B2,  the 
pitching  freedom  of  the  supporting  body  must  be  translated  into  that  of  the  single-cable 
attachment.  This  may  be  achieved  by  setting  the  inertia  of  the  supporting  body  to  zero 
while  holding  the  mass  constant.  If  this  mass  is  finite,  a  sliding-body  system  results,  in 
which  the  supporting  body  is  free  to  transverse  along  the  horizontal.  If  the  mass  is  infinite, 
the  system  simplifies  to  a  fixed-body  approximation. 


Sliding  Body 

Setting  Js  =  0  for  the  single-cable  system,  and  denoting  the  lengths  t  =  r  and  L  =  £+ a, 
the  characteristic  equation  becomes 

J/Aj4  -  ( migaL  +  Jigkm)u2  +  mig2 akm  =  0  (B16) 


with  roots 


U! 


2 


(migaL  +  Jigkm)  ±  [(migaL  +  Jigkm)2  -  4 Jiimig2akm}1/2 

2  Jil 


(B17) 
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Figure  B2:  Single-cable  Pendulum  System 


Fixed  Body 

Again,  setting  Js  —  0,  but  now  let  ms  — >  oo,  which  implies  that  km  — »  1, 

Jilco4  —  ( migaL  4-  Jig)ca 2  +  mig2a  =  0 
2  _  ( mtgaL  +  Jig)  ±  [{migaL  +  Jig)2  -  AMm^a}1/2 

u  -  w 

which  corresponds  to  the  solution  for  a  single  cable-suspended  body,  one  type  of  the 
“double-pendulum”  system. 

Multi-cable  System 

To  simulate  a  multi-cable  slung-load  system,  it  is  possible  to  either  use  Equations 
(B13)  and  (B15),  or  further  simplify  to  constrain  rotation  of  the  supporting  body,  shown 
in  Figure  B3.  This  may  be  implemented  by  taking  the  limit  as  the  inertia  approaches 
infinity.  Once  again,  the  corresponding  mass  can  be  adjusted  to  simulate  either  a  sliding 
or  fixed-body  system. 


(B18) 

(B19) 


Sliding  Body 

For  Js  — y  oo,  the  characteristic  equation  becomes 

(—  +  -^-fcm)w4  -  —kmuj2  =  0 
r  mir  r 

with  roots 

2  _  Jn  ml9akm  1 
U  \  ’  mia2  +  Jikm  J 


(B20) 

(B21) 
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Figure  B3:  Multi-cable  Pendulum  System 


Fixed  Body 


As  above,  Js  — >  oo.  Also,  letting  ms  oo  results  in  k„ 

+  -2L)W‘  _  =  0 

r  m/r  r 


1, 


u> 


{o, 

[  mtaz  +  Ji  ) 


(B22) 

(B23) 


Clearly,  this  corresponds  to  the  solution  for  a  single  pendulum  system  with  mass  mi 
and  inertia  Ji. 
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Appendix  C:  CH-47D  Chinook  Helicopter 
Aerodynamic  Model 


The  aerodynamic  model  used  for  the  CH-47D  was  based  on  a  set  of  linearised  stability 
and  control  derivatives  taken  from  Reference  [27].  The  derivatives,  as  well  as  associated 
trim  conditions,  were  computed  from  a  series  of  simulation  runs  at  various  airspeeds  using 
the  NASA  Ames  Research  Center  simulation  model. 

In  order  to  make  use  of  this  data,  the  derivatives  in  the  model  were  interpolated  between 
neighbouring  sets  with  respect  to  the  current  airspeed.  The  resultant  aerodynamic  forces 
and  moments  were  then  evaluated  as  a  linear  function  of  the  state  variables.  In  state-space 
form,  the  forces  can  be  written 


F a!  Trih  =  Ax  +  Bu 


(Cl) 


where  Fa  is  the  resultant  force,  x  is  the  vector  of  state  variables,  and  A  and  B  are  matrices 
comprising  the  stability  and  control  derivatives  respectively.  Expanding  Equation  (Cl) 


F Ax/mh 

xv 

xw 

xp 

xq 

XT  ' 

u 

' 

X5a 

XSr 

Xgc  ' 

F Ay!  ni>h 

Yu 

Yv 

Yw 

Yp 

Yq 

V 

Ysb 

YSa 

Ygr 

Yic 

F Az  /  mfi 

Zv 

%w 

Zp 

zq 

Zr 

w 

+ 

ZSb 

Z$a 

Zgr 

Zgc 

Max/ Ixx 

Lu 

Lu 

L<w 

Lp 

Lq 

Lp 

p 

Lsb 

Lfia 

Lgr 

Lgc 

MAy/Iyy 

Mu 

Mv 

Mw 

Mp 

Mq 

Mr 

q 

Msb 

Msa 

Mgr 

Mgc 

.  Maz/Izz  . 

.  Nu 

Nv 

Nw 

Np 

Nq 

Nr  _ 

r 

.  NSb 

NSa 

Ngr 

Ngc  . 

The  trim  conditions,  denoted  by  subscript  ”0”,  and  aerodynamic  derivatives,  in  column 
format,  for  each  of  the  tabulated  airspeeds  from  0.1  to  130  KTAS  are  summarised  in  the 
following  pages. 
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Run  #33  :  Straight  &  level  trim  ;  0.1  KTAS  ;  97.0  ft  ;  33000  lb 


uo 

Vo 

w0 

<t>  0 

do 

00 

^>0 

5 ao 

5ro 

<00 

0.1000 

0.0000 

0.0000 

-0.4453 

6.6120 

0.0000 

-0.0057 

0.2262 

-0.0032 

5.7555 

It 

V 

W 

P 

Q 

r 

56 

5r 

X 

-0.0200 

-0.0005 

0.0301 

0.0426 

2.7807 

-0.1590 

0.0570 

0.0000 

-0.0000 

0.9816 

Y 

-0.0003 

-0.1070 

0.0022 

-2.8362 

-0.0073 

-0.3305 

-0.0005 

1.0917 

0.0099 

0.0657 

Z 

0.0311 

0.0040 

-0.2983 

-0.1252 

-0.2647 

-0.0525 

0.0499 

-0.0001 

0.0000 

-8.4737 

L 

-0.0005 

-0.0108 

0.0004 

-1.2795 

0.0871 

-0.0903 

-0.0346 

0.4863 

-0.1263 

-0.0170 

M 

0.0111 

0.0000 

0.0006 

0.0295 

-1.0973 

-0.2686 

0.3282 

0.0000 

0.0000 

-0.0014 

N 

0.0007 

0.0006 

0.0000 

-0.0148 

-0.1338 

-0.0892 

0.0541 

0.0097 

0.1927 

-0.0006 

Run  #37  : 

Straight  &  level 

trim  ;  ! 

20.0  KTAS  ;  97.5  ft  ;  33000  lb 

u0 

^0 

Wo 

00 

do 

00 

5&0 

5ao 

5-ro 

5c0 

19.900 

0.0000 

1.9000 

-0.3858 

5.4756 

0.0000 

-1.4901 

0.2162 

0.2492 

5.6003 

u 

V 

w 

P 

Q 

T 

sb 

5Q 

Sf 

5C 

X 

-0.0128 

0.0001 

0.0268 

0.0277 

2.8998 

-0.1112 

0.0422 

0.0001 

-0.0001 

0.7238 

Y 

0.0021 

-0.0460 

0.0044 

-1.5122 

-0.1552 

-0.2018 

0.0632 

1.0867 

-0.0077 

0.0534 

Z 

-0.0405 

0.0044 

-0.2824 

-0.1057 

-2.7731 

0.1764 

0.2517 

-0.0010 

0.0012 

-8.4539 

L 

0.0004 

-0.0088 

0.0008 

-0.8862 

0.0239 

-0.0647 

-0.0105 

0.4856 

-0.1310 

-0.0218 

M 

0.0156 

-0.0009 

0.0166 

0.0171 

-1.6339 

-0.2477 

0.3355 

-0.0000 

0.0002 

-0.0083 

N 

0.0009 

0.0002 

0.0007 

-0.0059 

-0.1519 

-0.0764 

0.0525 

0.0086 

0.1917 

0.0021 
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Run  #41  :  Straight  k  level  trim  ;  40.0  KTAS  ;  98.0  ft  ;  33000  lb 


u0 

uo 

w0 

00 

00 

00 

h0 

ho 

ho 

5Co 

39.900 

0.0000 

3.0000 

-0.2941 

4.2599 

0.0000 

-2.3407 

0.2137 

0.3750 

5.1744 

u 

V 

W 

P 

q 

r 

h 

6a 

h 

X 

-0.0145 

0.0007 

0.0285 

0.0261 

2.8099 

-0.0915 

0.0353 

0.0000 

0.0000 

0.4914 

Y 

-0.0005 

-0.0594 

0.0071 

-1.7385 

-0.2622 

-0.2412 

0.1073 

1.0845 

-0.0233 

0.0654 

Z 

-0.1080 

0.0016 

-0.3498 

-0.0939 

-3.2915 

0.0318 

0.6498 

-0.0004 

0.0002 

-7.9678 

L 

-0.0003 

-0.0090 

0.0027 

-0.9513 

-0.0703 

-0.0839 

0.0163 

0.4856 

-0.1353 

-0.0096 

M 

0.0007 

-0.0015 

0.0259 

0.0128 

-1.7676 

-0.2562 

0.3744 

0.0000 

0.0001 

0.1235 

N 

0.0003 

0.0003 

0.0000 

-0.0048 

-0.0823 

-0.0683 

0.0370 

0.0076 

0.1913 

0.0054 

Run  #57a 

:  Straight  k  level  trim  ; 

60.0  KTAS  ;  98.6 

ft  ;  33000  lb 

u0 

^0 

w0 

00 

#0 

00 

O 

5(20 

ho 

ho 

60.000 

0.0000 

3.1000 

-0.2445 

2.9602 

0.0000 

-1.8241 

0.1971 

0.2589 

4.7454 

u 

V 

w 

P 

q 

r 

5fc 

5a 

5r 

5C 

X 

-0.0089 

0.0005 

0.0343 

0.0213 

2.6308 

-0.0700 

0.0475 

0.0000. 

0.0000 

0.3276 

Y 

-0.0014 

-0.0723 

0.0044 

-1.9788 

-0.1340 

-0.2488 

0.0941 

1.0852 

-0.0444 

0.0633 

Z 

-0.0755 

0.0024 

-0.5636 

-0.0627 

0.6404 

-0.1804 

0.8431 

-0.0004 

0.0003 

-8.1505 

L 

-0.0008 

-0.0093 

0.0020 

-1.0201 

-0.0664 

-0.0888 

0.0205 

0.4867 

-0.1417 

-0.0018 

M 

-0.0073 

-0.0013 

0.0145 

0.0199 

-1.5813 

-0.2736 

0.4118 

0.0000 

0.0001 

0.2338 

N 

0.0003 

0.0007 

-0.0003 

-0.0045 

-0.0501 

-0.0632 

0.0291 

0.0063 

0.1913 

0.0046 

Run  #5  7b 

:  Straight  k  level  trim  ; 

80.0  KTAS  ;  98.6 

ft  ;  33000  lb 

u0 

^0 

Wo 

00 

0o 

00 

h0 

5ao 

ho 

5c0 

80.000 

0.0000 

4.1000 

-0.2420 

2.9508 

0.0000 

-1.1228 

0.1754 

0.1152 

4.6284 

u 

V 

W 

P 

q 

r 

h 

5a 

h 

5c 

X 

-0.0057 

-0.0007 

0.0428 

0.0206 

2.6311 

-0.0647 

0.0499 

0.0000 

0.0000 

0.4185 

Y 

-0.0009 

-0.0878 

0.0024 

-2.0591 

-0.0412 

-0.2442 

0.0690 

1.0864 

-0.0310 

0.0473 

Z 

0.0230 

0.0049 

-0.6368 

-0.0530 

-0.4078 

-0.1833 

0.7232 

-0.0006 

0.0002 

-9.3412 

L 

-0.0007 

-0.0102 

0.0018 

-1.0358 

-0.0584 

-0.0864 

0.0182 

0.4865 

-0.1378 

-0.0032 

M 

-0.0081 

-0.0008 

0.0114 

0.0264 

-1.6518 

-0.2755 

0.4302 

0.0000 

0.0001 

0.2260 

N 

0.0003 

0.0013 

-0.0006 

-0.0118 

-0.0279 

-0.0622 

0.0251 

0.0071 

0.1915 

0.0050 
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Run  #57c  :  Straight  &  level  trim  ;  100.0  KTAS  ;  98.7  ft  ;  33000  lb 


Uo 

VO 

Wq 

4>  o 

#0 

6bo 

6ao 

6ro 

6co 

100.00 

0.0000 

4.6000 

-0.2679 

2.6354 

0.0000 

-0.6474 

0.1746 

0.0082 

4.7818 

u 

V 

W 

P 

Q 

r 

6b 

6a 

6r 

X 

-0.0126 

-0.0014 

0.0498 

0.0215 

2.6825 

-0.0649 

0.0500 

0.0001 

0.0000 

0.4893 

Y 

-0.0008 

-0.1051 

0.0008 

-2.0089 

0.0256 

-0.2401 

0.0566 

1.0891 

-0.0203 

0.0572 

Z 

0.0780 

0.0066 

-0.6769 

-0.0486 

-0.5284 

-0.1805 

0.6178 

-0.0008 

0.0003 

-10.3410 

L 

-0.0007 

-0.0113 

0.0013 

-1.0130 

-0.0210 

-0.0829 

0.0091 

0.4868 

-0.1350 

-0.0097 

M 

-0.0068 

-0.0004 

0.0117 

0.0321 

-1.7096 

-0.2839 

0.4471 

0.0000 

0.0002 

0.1998 

N 

0.0001 

0.0018 

-0.0007 

-0.0178 

-0.0608 

-0.0606 

0.0379 

0.0078 

0.1921 

0.0082 

Run  #57d 

:  Straight  &  level  trim  ; 

120.0  KTAS  ;  99.0  ft  ;  33000  lb 

uo 

vo 

Wo 

(f>0 

So 

TpO 

6 bo 

6ao 

6ro 

6 co 

120.10 

0.0000 

4.2000 

-0.3341 

1.9842 

0.0000 

-0.3143 

0.1934 

-0.1074 

5.1578 

u 

V 

W 

P 

q 

r 

6b 

<5j- 

6C 

X 

-0.0277 

-0.0016 

0.0523 

0.0237 

2.7440 

-0.0701 

0.0528 

-0.0000 

-0.0001 

0.5153 

Y 

-0.0001 

-0.1235 

0.0005 

-1.8659 

0.1258 

-0.2030 

0.0338 

1.0958 

-0.0102 

0.0643 

Z 

0.0631 

0.0082 

-0.7008 

-0.0400 

-0.4733 

-0.2141 

0.5106 

-0.0009 

0.0002 

-10.9880 

L 

-0.0005 

-0.0125 

0.0011 

-0.9612 

0.0029 

-0.0699 

0.0045 

0.4884 

-0.1328 

-0.0128 

M 

-0.0041 

0.0001 

0.0121 

0.0389 

-1.7035 

-0.2984 

0.4511 

0.0000 

0.0002 

0.1824 

N 

0.0001 

0.0021 

-0.0007 

-0.0256 

-0.0641 

-0.0669 

0.0412 

0.0085 

0.1933 

0.0135 

Run  #53  : 

:  Straight 

&  level 

trim  ; 

130.0  KTAS  ;  99.6 

ft  ;  33000  lb 

Uo 

Vo 

w0 

4>  0 

#0 

'‘I’O 

6bo 

6ao 

6ro 

6co 

130.20 

0.0000 

1.7000 

-0.3833 

0.7495 

0.0000 

-0.2349 

0.2098 

-0.1583 

5.4493 

U 

V 

w 

P 

q 

T 

6a 

<*c 

X 

-0.0378 

-0.0018 

0.0426 

0.0248 

2.7638 

-0.0789 

0.0658 

0.0000 

0.0000 

0.3553 

Y 

0.0003 

-0.1351 

0.0029 

-1.7771 

0.1674 

-0.1854 

0.0272 

1.1066 

-0.0152 

0.0593 

Z 

0.0132 

0.0082 

-0.7058 

-0.0358 

-0.3024 

-0.2564 

0.4628 

-0.0004 

0.0001 

-11.1430 

L 

-0.0005 

-0.0132 

0.0013 

-0.9303 

0.0149 

-0.0707 

0.0021 

0.4918 

-0.1355 

-0.0147 

M 

-0.0022 

0.0003 

0.0124 

0.0420 

-1.6772 

-0.3087 

0.4460 

0.0000 

0.0001 

0.1911 

N 

-0.0001 

0.0023 

-0.0007 

-0.0293 

-0.0689 

-0.0565 

0.0440 

0.0083 

0.1952 

0.0150 
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Appendix  D:  Matlab  Source  Code 


HSLJNIT.M 

7,  HSL_INIT  :  Initialisation  script  for  helicopter  slung-load  simulation 

7. 

7.  HSL.INIT 

7.  For  details  of  the  simulation,  see  reference  [1]  . 

7. 

7.  1 .  Stuckey ,  R .  A 

7.  "Mathematical  Modelling  of  Helicopter  Systems" 

7.  DST0-TR-000,  Defence  Science  and  Technology  Organisation,  Sep,  1999 

7.  R.A.  Stuckey  17/08/99  (c)  1999,  Defence  Science  and  Technology  Organisation 

7, - 

global  opt_  HDAT_  LDAT_  CDAT_  7.  GLOBAL  variables  (defined  below) 

7.  First,  define  the  simulation  parameters:  number  of  simulation  points,  time, 

7.  control  inputs,  plus  a  few  options  such  as  the  simulation  ares,  angular 
7.  representation,  the  cable  elasticity  and  nonlinear  solution  flags. 

7.  Also  define  the  number  of  bodies  in  the  system  (slung-loads  +  helicopter)  . 

7.  Number  of  points  for  simulation 
7,  Time  vector 

7,  Matrix  of  time  &  control  input  vectors 
7.  TD  =  [  t  db  da  dr  dc  ] 

7. 

7.  db  :  longitudinal  stick  (in) 

7.  da  :  lateral  stick  (in) 

7.  dr  :  pedal  (in) 

7.  dc  :  collective  (in) 

7.  Simulation  axes  [’longitudinal ’  I  ’ lateral'  I  ’ combined’] 

7.  Angular  representation  [’euler’ I ’quaternion’] 

7.  Elastic  cable  model  [Oil] 

7.  Nonlinear  simulation  [0 1 1] 

7.  Number  of  bodies 


N  =  101; 

t  =  [0:N-1]  ’*0.1; 

TD  =  [t,zeros(N,4)] ; 


opt_  =  {  ’com’ 
'eul  ’ 

0 

1  }; 

n  =  3; 


7.  Next,  create/load  the  data  structures  for  the  helicopter  and  loads.  These 
7.  will  contain  all  the  relevent  aerodynamic,  mass,  inertial  and  geometric 
7,  information.  Before  starting,  define  some  temporary  variables  (which  are  not 
7.  required  by  the  main  program)  for  convinience.  For  each  vector  here,  the 
7.  first  element  not  used  but  present  merely  to  maintain  consistency  in  the 
7,  numbering. 

eta  =  [  0  0.33  0.33  ]  ;  7.  Vector  of  mass  ratios 

lonL  =  [  0  0.0  0.6  ];  7.  Vector  of  length  ratios  (’O’  for  multiple-cable) 

L  =  [  0  20.0  25.0  ]  ;  7.  Vector  of  total  lengths  (ft) 

K  =  [  0  0.0  0.0  ];  7.  Vector  of  cable  stiffness  (’O’  for  inelasticity) 

7.  The  functions  ’hsl_ch47bdat ’  and  ’hsl.boxdat’  have  been  written  to  simplify 
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'/.  the  creation  of  standard  helicopter  and  load  data  structures.  The  inputs 
'/,  required  by  the  load  function  are,  respectively:  the  load  mass;  the  load 
’/,  inertia  matrix;  the  bridle  (sling-leg)  length;  the  sling  configuration;  and 
l  the  box  size  (ft)  .  Empty  matrices  given  for  any  of  these  will  trigger  their 
l  default  values.  The  data  structures  for  both  helicopter  and  loads  have 
'/,  exactly  the  same  fields. 


HDAT_ 

=  hsl_ch47bdat ;  7,  GLOBAL  struct  for  the  helicopter  and  load,  below 

LDAT_ 

=  struct  ( ’Name’  ,  []  ,  ... 

7.  Name  of  body 

’V’, 

7.  Velocity  vector  (kn) 

’A’,  [],... 

7.  Aerodynamic  derivative  matrix 

’m’,  □,  ... 

7,  Body  mass  (slugs) 

’J\  [],  ... 

7.  Inertia  matrix  (slugs. ft"2) 

’S’,  □); 

7,  Geometric  model  struct 

LDAT_ (2)  =  hsl_boxdat(eta(2)*HDAT_.m,  [],[], ’multiple’ ,  [  3.0  3.0  3.0  ]); 

LDAT_ (3)  =  hsl_boxdat(eta(3)*HDAT_.m, [] , (lonL(3)-l)*L(3) , ’single’, [  3.0  3.0  3.0  ]); 

'/,  The  geometric  model  struct  comprises  patches,  lines  and  links  for  rendering 
'/,  an  image  of  the  model.  The  links  are  also  used  in  the  mathematical  model  for 
7,  solution  of  the  equations  of  motion.  That  is,  the  (possibly  elastic)  slings 
7,  are  assumed  to  extend  from  the  link(s)  on  each  load  to  the  correspondign  link 
7,  or  attachment  point  on  the  helicopter.  These  structures  look  like: 

7. 

7,  S. Patches  .Data  :  Patch  matrix  (ft) 

7.  S. Patches .Coordlndex  :  Patch  coordinates 
7.  S. Lines. Data  :  Line  matrix  (ft) 

7.  S. Links. Data  :  Link  matrix  (ft) 

7,  Create  the  data  structures  for  the  cables  supporting  each  load.  These  contain 
7,  the  respective  link  attachment  indices,  cable  lengths,  cable  stretch  rates 
7,  and  the  stiffness  and  damping  constants.  The  link  attachment  matrix  lists  the 
7.  helicopter  link  indices  in  the  first  row  and  the  corresponding  load  link 
7.  indices  in  the  second  row  for  each  cable. 

CDAT_  =  struct  ( ’  i  ’  ,  []  ,  ...  7.  Link  attachment  matrix 


10’ 

. .  7.  Cable  lengths  - 

unloaded  (ft) 

1’ 

. .  7.  Cable  lengths  - 

current  (ft) 

ldot’ 

.  .  7.  Cable  stretch  rates  -  current  (ft/s) 

K’ 

,[]); 

7,  Cable  stiffness 

ft  damping  matrix 

7.  The  first  load  is  supported  by  multiple  (4)  cables  attached  to  the  first 
7.  (forward-most)  link  under  the  helicopter.  The  individual  cable  lengths  are 
7.  calculated  from  the  bridle  lengths  and  the  total  length  (helicopter 
7.  attachment  to  load  eg).  The  stiffness  matrix  is  empty,  since  this  simulation 
7,  is  inelastic  -  as  defined  in  opt_,  above  -  and  the  remaining  lengths  and 
7.  rates  are  set  to  their  initial  values. 

CDAT_(2)  .  i  =  [1111;1234]; 

CDAT_(2).10  =  sqrt (sum( (-LDAT_(2) .S. Links. Data+ones (4,1) * [  0  0  -L(2)  3). “2, 2))’; 
CDAT_(2).K  =  [];  CDAT_(2).l  =  CDAT_(2).10;  CDAT_  (2)  .  ldot  =  0*CDAT_(2)  .10; 

7,  The  second  load  is  supported  by  a  single  cable  attached  to  the  second  (middle) 
7,  helicopter  link.  The  cable  length  is  determined  by  the  cable-to-total  sling 
7.  length  ratio  and  the  other  fields  are  set  as  in  the  first  load. 

CDAT_  (3) .  i  =  [  3  ;  1  ]; 
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CDAT_(3) .10  =  lonL(3)*L(3) ; 

CDAT_(3).K  =  [];  CDAT_(3).l  =  CDAT_(3).10;  CDAT_(3)  .ldot  =  0*CDAT_(3)  .10; 


7,  The  next  step  is  to  compile  the  combined  system  mass/inertia  matrix  from  the 
7,  corresponding  helicopter  and  slung  load  matrices  defined  above .  To  constrain 
7.  the  motion  of  either  helicopter  or  load  in  any  of  the  [  x  y  z  ]  axes,  just 
7.  set  the  appropriate  element  (s)  of  the  explicitly  defined  [  0  0  0  ]  vectors 
7.  below  to  ’inf’.  Note  that  this  applies  to  both  rectillinear  motion,  for  the 
7.  mass  sub-matrices,  and  angular  motion  for  the  inertial  sub-matrices. 

D  =  zeros (n*6) ; 

jj  =  [1:3]  ; 

DC j j j)  =  HDAT_.m*eye(3)  +  diag([  0  0  0  ]);  7.  Helicopter  [  x  y  z  ]  mass 

D(n*3+jj  ,n*3+j  j)  =  HDAT_.  J  +  diag([  0  0  0  ]);  7.  Helicopter  [  x  y  z  ]  inertia 

for  j  =  2:n 

jj  =  (j-1)  *3+  [1:3]; 

DC  j  j ,  jj)  =  LDAT_(j)  .m*eye(3)  +  diag([  0  0  0  ]);  7.  Load  [  x  y  z  ]  mass 

D(n*3+j  j  ,n*3+j j)  =  LDAT_(j).J  +  diag([  0  0  0  ]);  7.  Load  [  x  y  z  ]  inertia 

end 


7.  Set  the  trim  state  for  the  helicopter  slung-load  system.  The  vector  uO 
7.  contains  the  generalised  velocity  coordinates.  The  vector  rO  contains  the 
7.  configuration  position  coordinates.  For  the  single  cable  case,  the  matrix 
7,  acO  contains  the  cable  angles  (in  cable-axes) .  These  variables  are  all 
7,  temporary  -  the  only  vector  required  by  the  main  program  is  the  trim  state 
7,  of  the  system  xO. 

uO  =  zeros(n*6,l) ;  7.  Trim  rate  vector  :  Single-cable  case 
7.  u  =  [  V1*_N  ;  Va22_c2  ;  .  .  .  ;  Vann_cn 
7.  wl_l  ;  w2_2  ;  . . .  ;  wn_n  ] 

7,  :  Multiple-cable  case 

7.  u  =  [  V1*_N  ;  Va22*_2  ;  .  .  .  ;  Vann*_n 

7.  wl_l  ;  w2_2  ;  ...  ;  wn_n  ] 

rO  =  zeros(n*6,l) ;  7.  Trim  position  vector 

7.  r  =  [  R1*_N  ;  R2*_N  ;  . . .  ;  Rn*_N 

7,  al  ;  a2  ;  . .  .  ;  an  ] 

acO  =  zeros (n, 3);  7.  Trim  cable-angle  matrix  :  Single-cable  case  only 

7.  ac  =  [  [  0  0  0  ]’  ac2  . . .  acn  ] 

7,  From  the  documentation  [*]  ,  the  sub-vectors  are  defined: 

7.  ... 

7,  V1*_N  :  helicopter  [x,y,z]  velocity  vector  in  inertial  axes 

7,  Vajj_cj  :  load  velocity  vectors  in  cable  axes 

7.  wl_l  :  helicopter  [p,q,r]  angular  velocity  vector  in  body  axes 
7.  wj_j  :  load  [p,q,r]  angular  velocity  vector  in  body  axes 
7,  Vajj*_j  :  load  velocity  vectors  in  load  axes 

7. 

7,  R1*_N  :  helicopter  [x,y,z]  position  vector  in  inertial  axes 

7.  Rj*_N  :  load  [x,y,z]  position  vector  in  inertial  axes 
7,  al  :  helicopter  [phi, theta, psi]  Euler  angular  position  vector 

7,  aj  :  load  [phi .theta, psi]  Euler  angular  position  vector 

7. 
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'/.  acj  :  cable/ cg-line  [phi, theta, psi]  angles  in  cable-axes 

’/,  The  function  ’hsl_load’  can  be  used  to  calculate  the  positions  of  the  loads, 
l  Rj*_N,  from  the  helicopter  positions,  the  helicopter  and  load  attitudes 
’/,  (all  included  in  the  vector,  rO) ,  the  cable  angles  (in  the  vector,  acO)  and 
I  the  cable  lengths  (supplied  globally  through  CDAT_.10). 

rO  =  hsl_load(rO,acO) ; 

'/.  Finally,  the  complete  state  vector  can  be  constructed  from  the  velocity  and 
'/.  position  vectors. 

xO  =  [uO;rO]  ; 


7,  Set  the  initial  state  for  the  helicopter  slung-load  system.  The  operations 
7.  are  exactly  the  same  as  those  for  the  trim  state  variables.  Setting  them  to 
7.  different  (non-zero)  values  will  create  an  initial  disturbance  in  the  system. 

u  =  zeros (n*6,l) ; 
u(l:3)  =  [  0  0  0  ]'; 


r  =  zeros(n*6,l) ; 
r(n*3+[l:n*3] )  =  [  0  0 

-5  10 

0  0 

ac  =  acO; 

ac(:,l:n)  =  [  0  0 

0  0 
5  -10 


0  ... 

0  ... 

0  ] ’*pi/180; 

0 

0 

0  3 ’*pi/180; 


r  =  hsl_load(r,ac) ; 


x  =  [u;r]  ; 


7,  Now  clear  extraneous  variables  from  the  workspace  and  run  the  main  program. 
7.  The  required  inputs,  as  specified  above,  axe: 

7. 

7.  opt_  HDAT_  LDAT_  CDAT_  n  t  TD  D  xO  x 
7. 

7.  The  outputs  are: 

7. 

7,  U  :  Generalised  velocity  coordinate  matrix 

7,  R  :  Configuration  position  coordinate  matrix  using  Euler  angles 

7.  R_q  :  Configuration  position  coordinate  matrix  using  quaternions 

7.  X  :  Full  state  variable  matrix  using  Euler  angles 

7.  X_q  :  Full  state  variable  matrix  using  quaternions 

clear  eta  lonL  L  K  j  jj  uO  rO  acO  u  r  ac 

hslsim 

7,  Lastly,  Plot  the  primary  state  and  control  variables  as  time  histories, 
xn  =  [i:N];  hslplot 
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HSL_CH47BDAT.M 


function  D  =  hsl_ch47bdat (mass .inertia, links) 

'/.  HSL_CH47BDAT  :  Aircraft  aerodynamic,  mass,  inertial  and  geometric  data 

7. 

7.  DAT  =  HSL_CH47BDAT((mass)  (.inertia)  (.links)) 

7. 

7.  mass  :  Helicopter  mass  (slugs) 

7,  inertia  :  Inertia  matrix  (slug.ft"2) 

7,  links  :  Links  data  matrix  (see  HSL_INIT) 

7. 

7,  DAT  :  Data  structure  (see  HSL.INIT) 

7.  Hass,  inertial  and  aerodynamic  data  obtained  from  reference  [1]. 

7. 

7,  1.  Weber,  J.M.  and  Liu,  T.Y.  and  Chung,  W. 

7.  "A  Mathematical  Simulation  Model  of  the  CH-47B  Helicopter,  Volumes  1  &  2" 

7.  NASA-TM-84351-V0L-1/2,  NASA  Ames  Research  Center,  Aug,  1984 

’l  R.A.  Stuckey  01/06/99  (c)  1999,  Defence  Science  and  Technology  Organisation 

7. - 


if  nargin<3,  links  =  []  ;  end 

if  nargin<2,  inertia  =  []  ;  end 

if  nargin<l,  mass  =  []  ;  end 

7,  Set  some  defaults 

if  isempty (mass) ,  mass  =  33000/32.174;  end 

if  isempty (inertia) 

inertia  =  [  34000  0  14900 

0  202500  0 

14900  0  191000  ]; 

end 

D  =  struct( ’Name’ ,  □  , ’V’ ,  □  ,  ’  A’ ,  []  , ’m’ ,  □  , '  J’ ,  []  , ’S  ’ ,  []  )  ; 


D .Name  =  ’CH-47B’  ; 

D.V  =  [  0.1  20  40  60  80  100  120  130  ]’;  7.  Trim  airspeed  (KTAS) 


7,  Construct  the  aerodynamic  coefficient  matrix 
D.A  =  zeros (8, 10*6,2) ; 


7.  SAS  Off 
7. 

7.  v  =  0.1  KTAS 

7. 

A  =  [ 

-1.99980e-002 
-4 . 59450e-004 
3 . 00850e-002 
4.26000e-002 
2.78070e+000 
-1.59040e-001 
5.70000e-002 


-2.78070e-004 
-1.07040e-001 
2 . 23380e-003 
-2.83620e+000 
-7.31250e-003 
-3.30450e-001 
-4.84860e-004 


3. 11490e-002 
4.00190e-003 
-2.98310e-001 
-1 . 25200e-001 
-2.64720e-001 
-5.25060e-002 
4.98650e-002 


-5.32570e-004 
-1 . 07710e-002 
3.99760e-004 
-1 . 27950e+000 
8.71280e-002 
-9.02880e-002 
-3.46400e-002 


1 . 10900e-002 
9 . 49780e-006 
6 . 38150e-004 
2.95300e-002 
-1.09730e+000 
-2.68610e-001 
3 . 28160e-001 


7.34520e-004 

5.81320e-004 

4.24930e-005 

-1.48080e-002 

-1.33780e-001 

-8.91680e-002 

S.40720e-002 
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1.00060e-005 
-4.76460e-007 
9.81570e-001 
]’;  D.A(1, :  ,1)  < 

7.  v  =  20.0  KTAS 
7. 

A  =  [ 

-1.28440e-002 
1.43630e-004 
2.67660e-002 
2.76530e-002 
2.89980e+000 
-1.11160e-001 
4.21820e-002 
9.33870e-005 
-7.28990e-005 
7.23840e-001 
]’;  D.A(2, :  ,1)  = 

7.  v  =  40.0  KTAS 
7. 

A  =  [ 

-1.45080e-002 
7 . 12390e-004  • 
2.85410e-002 
2 . 60640e-002  • 
2 . 80990e+000 
-9 . 15160e-002 
3 . 53470e-002 
3 . 09700e-005 
9 . 52930e-006  ■ 
4 . 91350e-001 
]>;  D .  A(3, : ,  1)  = 

7.  v  =  60.0  KTAS 
7. 

A  =  [ 

-8.90430e-003  • 
4.56690e-004  ■ 
3.42780e-002 
2 . 12800e-002  ■ 
2 . 63080e+000  ■ 
-6 . 99840e-002  • 
4.75450e-002 
2.38230e-005 
9.52930e-006  ■ 
3.27550e-001 
]>;  D.A(4, :  ,1)  = 

7.  v  =  80.0  KTAS 
7. 

A  =  [ 

-5.70700e-003  • 
-7.31440e-004  • 
4 . 27800e-002 
2 . 05520e-002  • 
2.63110e+000  ■ 


1 . 09170e+000 
9. 86840e-003 
6. 56800e-002 
;  A( : )  ’ ; 


2.05350e-003 
•4.60310e-002 
4.35330e-003 
■1 . 51220e+000 
•1 . 55240e-001 
•2.01820e-001 
6.31740e-002 
1.08670e+000 
7.71500e-003 
5.33870e-002 
’  A  ( : )  ’ ; 


4. 56290e-004 
5.94030e-002 
7.10780e-003 
■1.73850e+000 
2. 62180e-001 
2.41250e-001 
1.07320e-001 
1.08450e+000 
2.32600e-002 
6 . 53500e-002 
;  A( : )  ’ ; 


1.36390e-003 
7 . 23000e-002 
4.36120e-003 
1 . 97880e+000 
1.34030e-001 
2.48750e-001 
9.41300e-002 
1.08520e+000 
4.43990e-002 
6.33380e-002 
:  A( : )  ’  ; 


8.59470e~004 
8.77940e-002 
2.40690e-003 
2.05910e+000 
4. 12240e-002 


-7.62340e-005 
0 . 00000e+000 
-8 . 47370e+000 


-4.05230e-002 
4 . 43150e-003 
-2 . 82440e-001 
-1.05700e-001 
-2 . 77310e+000 
1.76410e-001 
2.51680e-001 
-1.04440e-003 
1.22740e-003 
-8.45390e+000 


-1.07970e-001 
1 . 55100e-003 
-3.49780e-001 
-9 . 39210e-002 
-3.29150e+000 
3. 18470e-002 
6 . 49750e-001 
-4 . 11670e-004 
1.90590e-004 
-7.96780e+000 


-7 . 54850e-002 
2 . 38840e-003 
-5 . 63570e-001 
-6 . 27030e-002 
6 . 40370e-001 
-1 . 80370e-001 
8 . 43060e-001 
-3.58300e-004 
2 . 59200e-004 
-8 . 15050e+000 


2 . 30150e-002 
4 . 89230e-003 
-6 . 36790e-001 
-5 . 30210e-002 
-4 . 07800e-001 


4.86300e-001 
-1 . 26340e-001 
-1 . 70400e-002 


3.66750e-004 
-8.76990e-003 
8.33700e-004 
-8.86220e-001 
2.38660e-002 
-6.46980e-002 
-1.04610e-002 
4.85580e-001 
-1.30970e-001 
-2. 18100e-002 


-2 . 78280e-004 
-8.99530e-003 
2.74400e-003 
-9.51320e-001 
-7.03350e-002 
-8.39020e-002 
1 . 62620e-002 
4.85580e-001 
-1.35310e-001 
-9.56590e-003 


-7.78620e-004 
-9.29200e-003 
2.04470e-003 
-1 . 02010e+000 
-6 . 64380e-002 
-8.87950e-002 
2.05270e-002 
4.86690e-001 
-1.41650e-001 
-1.81260e-003 


-7.20310e-004 
-1.01780e-002 
1.81170e-003 
-1.03580e+000 
-5 . 84190e-002 


7.42500e-006 

0.00000e+000 

-1.43180e-003 


1.55770e-002 
-8.73300e-004 
1.66480e-002 
1.71360e-002 
-1.63390e+000 
-2.47710e-001 
3.35480e-001 
-2.16560e-005 
2. 10380e-004 
-8.30050e-003 


7.31330e-004 
-1.53180e-003 
2.58770e-002 
1 . 28110e-002 
-1.76760e+000 
-2.56190e-001 
3.74440e-001 
7.42500e-006 
1.06420e-004 
1.23530e-001 


-7.26000e-003 
-1.28440e-003 
1.45080e-002 
1.99240e-002 
-1 . 58130e+000 
-2.73610e-001 
4 . 11840e-001 
1.48500e-005 
1.04570e-004 
2.33770e-001 


-8.07210e-003 
-7.76930e-004 
1 . 14460e-002 
2.64380e-002 
-1.65180e+000 


9 . 74910e-003 
1 . 92690e-001 
-5.52630e-004 


9.33590e-004 
2 . 30440e-004 
7 . 13630e-004 
-5.86890e-003 
-1 . 51860e-001 
-7 . 64320e-002 
5.25400e-002 
8 . 61240e-003 
1.91730e-001 
2.06350e-003 


2.52740e-004 
3 . 20370e-004 
2.49070e-005 
-4.82000e-003 
-8 . 23200e-002 
-6 . 83260e-002 
3.69750e-002 
7 . 62330e-003 
1 . 91260e-001 
5 . 36610e-003 


2.75770e-004 
7.42330e-004 
-3.09370e-004 
-4 . 54600e-003 
-5 . 01390e-002 
-6.32470e-002 
2.91390e-002 
6.29670e-003 
1 . 91280e-001 
4.56090e-003 


2.93460e-004 
1.34520e-003 
-5 . 56030e-004 
-1 . 17890e-002 
-2 . 79330e-002 
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-6.46630e-002  -2 . 44250e-001  - 
4.99400e-002  6.90110e-002 

4.24050e-005  1.08640e+000  - 

1 . 14350e-005  -3. 10170e-002 
4 . 18550e-001  4.72800e-002  - 

]>;  D.A(5, : ,1)  =  A(:)>; 

7.  v  =  100.0  KTAS 

7. 

A  =  [ 

-1.26360e-002  -7.51370e-004 
-1.35680e-003  -1.05140e-001 
4.98200e-002  8.41670e-004  - 

2.15280e-002  -2.00890e+000  - 
2.68250e+000  2.56190e-002  - 

-6.49300e-002  -2.40140e-001  - 
5.00190e-002  5.66120e-002 

6.24170e-005  1.08910e+000  - 

1 . 04820e-005  -2 . 02930e-002 
4.89300e-001  5.72410e-002  - 

]’;  D.A(6,:,1)  =  A(:)’; 

7.  v  =  120.0  KTAS 
7. 

A  =  [ 

-2.77000e-002  -1.34170e-004 
-1.64920e-003  -1.23550e-001 
5.22720e-002  5.28050e-004  - 

2.36600e-002  -1.86590e+000  - 
2.74400e+000  1.25820e-001  - 

-7.00630e-002  -2.03000e-001  - 
5.27850e-002  3.37610e-002 

-2.66820e-005  1.09580e+000  - 

-7.33760e-005  -1.02440e-002 
5 . 15320e-001  6.42760e-002  - 

]>;  D.A(7 , : ,1)  =  A(:)>; 

7.  v  =  130.0  KTAS 
7. 

A  =  [ 

-3.78490e-002  2.51630e-004 

-1.78130e-003  -1.35140e-001 
4.26380e-002  2.86170e-003  - 

2.47760e-002  -1 .77710e+000  - 
2.76380e+000  1.67400e-001  - 

-7.88740e-002  -1 .85370e-001  - 
6.57840e-002  2.71970e-002 

2.38230e-005  1.10660e+000  - 

1.76290e-005  -1 . 51530e-002 
3.5S270e-001  5.92710e-002  - 

]’;  D .  A(8, : ,  1)  =  A(:)'; 


7.  SAS  On 
7. 

7.  v  =  0.1  KTAS 
7. 

A  =  [ 


1.83310e-001 
7. 23190e-001 
5.71760e-004 
2. 28700e-004 
9.34120e+000 


7.79880e-002 
6.60610e-003 
6.76910e-001 
4.85800e-002 
5 . 28360e-001 
1.80540e-001 
6. 17770e-001 
7.85210e-004 
3. 12560e-004 
1.03410e+001 


6 . 30760e-002 
8 . 15140e-003 
7. 00810e-001 
3.99660e-002 
4.73260e-001 
2 . 14120e-001 
5 . 10630e-001 
8 . 69070e-004 
1 . 67720e-004 
1 . 09880e+001 


1 . 32080e-002 
8. 17960e-003 
7. 05780e-001 
3. 57920e-002 
3. 02350e-001 
2. 56380e-001 
4. 62770e-001 
3 . 58300e-004 
5 . 33640e-005 
1 . 11430e+001 


-8.63580e-002 
1.81860e-002 
4.86460e-001 
-1.37820e-001 
-3 . 23580e-003 


-6.59260e-004 
-1.12800e-002 
1.32360e-003 
-1.01300e+000 
-2 . 10030e-002 
-8 . 29250e-002 
9.06130e-003 
4.86820e-001 
-1.34970e-001 
-9.66370e-003 


-4.79240e-004 
-1 . 24840e-002 
1 . 12150e-003 
-9.61180e-001 
2.85120e-003 
-6 . 99070e-002 
4.49670e-003 
4 . 88370e-001 
-1.32770e-001 
-1.27980e-002 


-4.70110e-004 
-1.31600e-002 
1.31020e-003 
-9.30250e-001 
1.48530e-002 
-7.06540e-002 
2 . 11230e-003 
4.91790e-001 
-1.35500e-001 
-1.47050e-002 


-2.75520e-001 

4.30170e-001 

1.98000e-005 

1.31180e-004 

2.26040e-001 


-6.81070e-003 
-3.57580e-004 
1 . 16510e-002 
3.20850e-002 
-1.70960e+000 
-2.83950e-001 
4.47050e-001 
1.36120e-005 
1.52830e-004 
1.99790e-001 


-4. 14310e-003 
9 . 20080e-005 
1 . 20600e-002 
3.89050e-002 
-1.70350e+000 
-2.98380e-001 
4.51130e-001 
8.04370e-006 
1.72630e-004 
1.82350e-001 


-2.21120e-003 
2.78040e-004 
1.24230e-002 
4. 19650e-002 
-1.67720e+000 
-3.08690e-001 
4.46010e-001 
6.80620e-006 
7 . 92000e-005 
1.91100e-001 


-6.21810e-002 
2.50860e-002 
7. 14560e-003 
1.91540e-001 
5.02230e-003 


8.58610e-005 

1.77010e-003 

-7.02930e-004 

-1.78010e-002 

-6.07910e-002 

-6.06080e-002 

3.78660e-002 

7.83530e-003 

1.92080e-001 

8.24070e-003 


8.77130e-005 
2. 10200e-003 
-7.09630e-004 
-2.55800e-002 
-6.41210e-002 
-6.68580e-002 
4. 12170e-002 
8.50380e-003 
1.93300e-001 
1.35210e-002 


-1.31000e-004 

2.27430e-003 

-6.89930e-004 

-2.93120e-002 

-6.89040e-002 

-5.64850e-002 

4.39730e-002 

8.25240e-003 

1.95170e-001 

1.50040e-002 
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-2.00190e-002 
-4.58720e-004 
2 -95890e-002 
4.26580e-002 
2.76530e+000 
-1.59080e-001 
5.69980e-002 
-4.76460e-007 
0 . OOOOOe+OOO 
9.82140e-001 
]';  D.AC1, :  ,2) 


-2 . 79300e-004 
-1.07040e-001 
2. 20110e-003 
-4.80530e+000 
-8.48820e-003 
-3 . 54870e-001 
-4 . 73250e-004 
1.09170e+000 
9.87340e-003 
6.57510e-002 
=  A(:>>; 


3. 13430e-002 
3.99430e-003 
-2 . 93800e-001 
-1 . 25310e-001 
-1 . 16850e-001 
-5 . 24300e-002 
4 . 98950e-002 
0. OOOOOe+OOO 
-7.62340e-006 
-8.47890e+000 


-4.05190e-002 
4.44410e-003 
-2.82440e-001 
-1.00800e-001 
-2.77560e+000 
1.74270e-001 
2 . 51470e-001 
2.82070e-004 
1 . 15880e-003 
-8.45410e+000 


-1.04060e-001 
1 . 55750e-003 
-3.49160e-001 
-9.25870e-002 
-3 . 29140e+000 
3.20760e-002 
6.49480e-001 
-4 . 11670e-004 
1 . 90590e-004 
-7 . 96770e+000 


-6 . 34290e-002 
2.40440e-003 
-5 . 62920e-001 
-6.09110e-002 
-6.40410e-001 
-1.80520e-001 
8.43140e-001 
-3 . 88800e-004 
1 . 82960e-004 
-8 . 15030e+000 


-5.33410e-004 
-1.07860e-002 
3.92040e-004 
-2. 16340e+000 
8.65030e-002 
5.98430e-002 
-3.46380e-002 
4.86300e-001 
-1.26340e-001 
-i.70130e-002 


3.68320e-004 
-1.19320e-002 
8.35320e-004 
-1.76990e+000 
2.38450e-002 
9.46100e-002 
-1 . 04580e-002 
4.85580e-001 
-1.30970e-001 
-2 . 18080e-002 


-1.66880e-004 
-1 . 53840e-002 
2.75790e-003 
-1.83520e+000 
-7.03380e-002 
-8.39040e-002 
1.62660e-002 
4.85580e-001 
-1 . 35310e-001 
-9 . 56410e-003 


-4.90270e-004 

-1.89950e-002 

2.06160e-003 

-1.90630e+000 

-6.64230e-002 

-8.87900e-002 

2.05240e-002 

4.86690e-001 

-1.41660e-001 

-1.81280e-003 


1.10900e-002 
8 . 69340e-006 
6 . 48540e-004 
2 . 95170e-002 
-1 . 09920e+000 
-2 . 68650e-001 
3. 28170e-001 
6. 18750e-007 
6. 18750e-007 
-1.43490e-003 


l,55760e-002 
-8.70520e-004 
1.66470e-002 
1.76080e-002 
-1.63400e+000 
-2.47990e-001 
3.35490e-001 
5 . 07370e-005 
1 . 16320e-004 
-8.30980e-003 


3.02760e-003 
-1 . 52910e-003 
2.61800e-002 
1 . 27460e-002 
-1 . 76760e+000 
-2 . 56190e-001 
3.74470e-001 
1.67060e-005 
1 . 01470e-004 
1 . 23540e-001 


-1 . 49840e-003 
-1 . 27800e-003 
1 . 48370e-002 
2.00350e~002 
-1.58130e+000 
-2.73610e-001 
4. 11830e-001 
4.33130e-006 
1 . 03330e-004 
2.33750e-001 


7.34570e-004 
6 . 02580e-004 
4 . 01240e-005 
-2.27710e-002 
-1.33820e-001 
-3.23810e-001 
5.40740e-002 
9.74810e-003 
1.92690e-001 
-5.47820e-004 


9.31240e-004 
4.54010e-003 
7 . 11290e-004 
-1 . 10830e-002 
-1.51840e-001 
-3 . 10110e-001 
5.25420e-002 
8 . 61630e-003 
1.91730e-001 
2 . 06790e-003 


4 . 71550e-004 
8 . 83450e-003 
5 . 38230e-005 
-8 . 27830e-003 
-8.23200e-002 
-6.83270e-002 
3.69740e-002 
7.62280e-003 
1.91260e-001 
5 . 36710e-003 


6 . 84190e-004 
1.33160e-002 
-2 . 86400e-004 
-5 . 59910e-003 
-5.01410e-002 
-6.32470e-002 
2.91370e-002 
6 . 29400e-003 
1.91280e-001 
4.55900e-003 


•/.  v 

7. 

A  = 
-1. 
1. 

2. 

2. 

2. 

-1. 

4. 

-1. 

-8. 

7. 


20.0  KTAS 


[ 

28450e-002 
43680e-004 
67660e-002 
73040e-002 
90000e+000 
11040e-001 
22050e-002 
95350e-005 
67170e-005 
23850e-001 
D.A(2, :  ,2) 


2.05360e-003 
-4.69260e-002 
4.35350e-003 
-3.47460e+000 
-1.55300e-001 
-1.93430e-001 
6.31900e-002 
1 . 08670e+000 
-7 . 72250e-003 
5.34030e-002 
=  A( : )  ’ ; 


7.  v 

7. 

A  = 
-1. 
7. 
2. 
2. 
2. 
-9. 
3. 

3. 
7. 

4. 


=40.0  KTAS 

[ 

42800e-002 
12670e-004 
85590e-002 
59480e-002 
80990e+000 
15410e-002 
53760e-002 
38290e-005 
62340e-006 
91340e-001 
D.A(3, : ,2) 


2.09760e-004 
-6 . 16240e-002 
7. 19340e-003 
-3.69780e+000 
-2.62190e-001 
-2.41250e-001 
1.07330e-001 
1 . 08450e+000 
-2.32630e-002 
6 . 53570e-002 
=  A(:)>; 


7.  v  =  60.0  KTAS 
7. 

A  =  [ 

-8.27460e-003 
4.55900e-004 
3.43190e-002 
2. 11980e-002 
2.63080e+000 
-6.99690e-002 
4.75450e-002 
1.90590e-005 
1.28650e-005 
3.27530e-001 
]’;  D.  A(4,  :  ,2) 


-4. 85040e-005 
-7.65220e-002 
4. 43620e-003 
-3.94060e+000 
-1.33990e-001 
-2.48740e-001 
9.41200e-002 
1.08520e+000 
-4.44340e-002 
6.33311e-002 
=  A(:)>; 


7.  v  =  80.0  KTAS 
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7. 

A  =  [ 

-5.01240e-003  1.05650e-004 

-7.36350e-004  -8.98880e-002 
4 . 28110e-002  2.45010e-003 

2.0S250e-002  -4.02230e+000 
2.63120e+000  -4. 12660e-002 
-6.46790e-002  -2.44270e-001 
4.99380e-002  6.90080e-002 

4 . 19290e-005  1.08640e+000 

1.71530e-005  -3.09520e-002 
4. 18540e-001  4.72930e-002 

]’;  D.  A(5, :  ,2)  =  A(:)'; 

7.  v  =  100.0  KTAS 
7. 

A  =  [ 

-1.21400e-002  -2. 12260e-004 
-1.35350e-003  -1 . 06330e-001 
4.98460e-002  8.65710e-004 

2.15560e-002  -3.97640e+000 
2.68240e+000  2.56070e-002 

-6.49100e-002  -2 ,40130e-001 
5.00630e-002  5.66060e-002 

5.47930e-005  1.08910e+000 

1.66760e-005  -2.03130e-002 
4.89350e-001  5.72900e-002 

]’;  D .  A(6, :  ,2)  =  A(:)’; 

7  v  =  120.0  KTAS 

7. 

A  =  [ 

-2.72310e-002  1.60600e-004 

-1.64710e-003  -1 . 24180e-001 
5 . 22970e-002  5.26710e-004 

2.33190e-002  -3.84480e+000 
2.74360e+000  1.11460e-001 

-6.97970e-002  -2.02890e-001 
5.26670e-002  3.35910e-002 

6.76580e-005  1.09570e+000 

2 . 14410e-005  -1 . 01920e-002 
5 . 15400e-001  6.40920e-002 

]>;  D.A(7, :  ,2)  =  A(:)>; 

7.  v  =  130.0  KTAS 
7. 

A  =  [ 

-3.72690e-002  4.92160e-004 

-1.78010e-003  -1 .35900e-001 
4.26400e-002  2.85720e-003 

2.47090e-002  -3.77600e+000 
2 . 76380e+000  1 . 67250e-00 1 

-7.88690e-002  -1 . 85400e-001 
6.57720e-002  2.72330e-002 

3.85940e-005  1.10660e+000 

1.76290e-005  -1.51090e-002 
3 . 55260e-001  5.91680e-002 

]’;  D.A(8,:,2)  =  A(:)’; 


3.32200e-002 

4.92400e-003 

-6.36340e-001 

-5.29830e-002 

-4.07570e-001 

-1.83170e-001 

7.23020e-001 

-5.94630e-004 

1.82960e-004 

-9.34120e+000 


8.45320e-002 
6.61450e-003 
-6.76670e-001 
-4.71510e-002 
-5.22640e-001 
-1 . 80540e-001 
6. 18090e-001 
-7.16600e-004 
2.97310e-004 
-1.03410e+001 


6.76800e-002 
7.99580e-003 
-7.00690e-001 
-3.88800e-002 
-4.81120e-001 
-2.21880e-001 
5 . 13600e-001 
-8.30950e-004 
1.52470e-004 
-1.09850e+001 


1.74470e-002 
8.34420e-003 
-7.05740e-001 
-3.S1820e-002 
-3.02540e-001 
-2.56340e-001 
4.62820e-001 
-5.26020e-004 
2.28700e-005 
-1 . 11430e+001 


-4.69290e-004 
-1.68720e-002 
1.82300e-003 
-1 . 92140e+000 
-5.84350e-002 
-8.63610e-002 
1.81850e-002 
4.86460e-001 
-1.37800e-001 
-3. 23290e-003 


-5. 70060e-004 
-1 . 68210e-002 
1 . 32940e-003 
-1.89920e+000 
-2. 10470e-002 
-8 . 29660e-002 
9.06090e-003 
4. 86810e-001 
-1.34970e-001 
-9. 64680e-003 


-4.40090e-004 
-1.74270e-002 
1. 12160e-003 
-1 . 84990e+000 
2.34480e-003 
-6. 98590e-002 
4.45150e-003 
4.88340e-001 
-1 . 32750e-001 
-1 . 28180e-002 


-4. 51350e-004 
-1 . 80350e-002 
1 . 30860e-003 
-1.82540e+000 
1.47950e-002 
-7.06680e-002 
2. 12550e-003 
4. 91790e-001 
-1 . 35490e-001 
-1.47390e-002 


-2.04460e-003 
-7.82160e-004 
1. 17170e-002 
2.64020e-002 
-1.65180e+000 
-2.75610e-001 
4.30170e-001 
8.04370e-006 
1.31790e-004 
2.26060e-001 


-2.20040e-003 

-3.52720e-004 

1.18410e-012 

3.21750e-002 

-1.70950e+000 

-2.84000e-001 

4.47080e-001 

1.67060e-005 

1.77580e-004 

1.99830e-001 


-1.53330e-004 
1.00270e-004 
1.21910e-002 
3.86730e-002 
-1.70360e+000 
-2.98370e-001 
4.51180e-001 
9.28120e-006 
1.75720e-004 
1 . 82370e-001 


1.73700e-003 
2 . 81930e-004 
1.24690e-002 
4.20970e-002 
-1.67700e+000 
-3.08690e-001 
4.46080e-001 
1 . 29940e-005 
8.78620e-005 
1.91010e-001 


6.48800e-004 
1.03840e-002 
-5.37880e-004 
-1.43840e-002 
-2.79310e-002 
-6.21850e-002 
2 . 50850e-002 
7. 14940e-003 
1.91540e-001 
5.02440e-003 


4.79490e-004 
9.49190e-003 
-6.89160e-004 
-2. 15830e-002 
-6.07840e-002 
-6.06040e-002 
3.78660e-002 
7.83420e-003 
1 . 92080e-001 
8.24110e-003 


4.51540e-004 
9. 18990e-003 
-7.00720e-004 
-3.06020e-002 
-6.39660e-002 
-6.68540e-002 
4. 11780e-002 
8.50630e-003 
1.93290e-001 
1.34750e-002 


2.56790e-004 
9.20090e-003 
-6.85800e-004 
-3.36860e-002 
-6.88860e-002 
-5.64830e-002 
4.39700e-002 
8 . 25660e-003 
1.95170e-001 
1.50090e-002 
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D.m  =  mass;  '/,  Mass  (slugs) 

D.J  =  inertia;  */.  Inertia  matrix  (slugs. ft"2) 


'/,  Construct  the  geometric  struct 

D.S  =  struct  ( ’Patches’ ,[], ’Lines’ ,[], ’Links’ ,[]  ) 


D.S. Patches. 

Data  =  [ 

27.57681 

1.60769 

4.56059 

23.11464 

4.75745 

6.95572 

23.11464 

4.75745 

-2.49356 

-25.42775 

4.75745 

-2.49356 

-13.96066 

6.69324 

6.89010 

-13.96066 

-6.69324 

6.89010 

-25.42775 

-4.75745 

-2.49356 

23.11464 

-4.75745 

-2.49356 

23.11464 

-4.75745 

6.95572 

27.57681 

-1.60769 

4.56059 

23.11464 

1.27959 

-2.46075 

23.11464 

1.27959 

-4.79026 

15.73240 

1.27959 

-4.42935 

15.73240 

1.27959 

-2.46075 

15.73240 

-1.27959 

-2.46075 

15.73240 

-1.27959 

-4.42935 

23.11464 

-1.27959 

-4.79026 

23.11464 

-1.27959 

-2.46075 

-8.21890 

1.27959 

-2.46075 

-10.41718 

1.27959 

-8.30093 

-25.08325 

1.27959 

-9.58052 

-25.42775 

1 . 27959 

-2.46075 

-25.42775 

-1.27959 

-2.46075 

-25.08325 

-1.27959 

-9.58052 

-10.41718 

-1.27959 

-8.30093 

-8.21890 

]; 

-1.27959 

-2.46075 

D. S .Patches . 

Coordlndex 

=  [ 

0  1 

2 

2  -1  . 

1  2 

3 

4  -1  . 

3  4 

5 

6  -1  . 

5  6 

7 

8  -1  . 

7  8 

9 

9  -1  . 

1  4 

5 

8  -1  . 

0  1 

8 

9  -1  . 

0  2 

7 

9  -1  . 

2  3 

6 

7  -1  . 

10 

11 

12 

13 

-1 

12 

13 

14 

15 

-1 

14 

15 

16 

17 

-1 

10 

11 

16 

17 

-1 

11 

12 

15 

16 

-1 

18 

19 

20 

21 

-1 

20 

21 

22 

23 

-1 

22 

23 

24 

25 

-1 
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18  19  24  25  -1  .  . . 

19  20  23  24  -1 

]; 


D.S. Lines. Data  =  [ 


51.47889 

0.00000 

-7.64473 

47.46121 

14.99417 

-7.64473 

36.48472 

25.97066 

-7.64473 

21.49055 

29.98834 

-7.64473 

6.49638 

25.97066 

-7.64473 

-4.48011 

14.99417 

-7.64473 

-8.49779 

0.00000 

-7.64473 

-4.48011 

-14.99417 

-7.64473 

6.49638 

-25.97066 

-7.64473 

21.49055 

-29.98834 

-7.64473 

36.48472 

-25.97066 

-7.64473 

47.46121 

-14.99417 

-7 . 64473 

51.47889 

-0.00000 

-7.64473 

NaN 

NaN 

NaN 

10.30234 

0.00000 

-13.22243 

6.28466 

14.99417 

-13.22243 

-4.69183 

25.97066 

-13.22243 

-19.68600 

29.98834 

-13.22243 

-34.68017 

25.97066 

-13.22243 

-45.65666 

14.99417 

-13.22243 

-49.67434 

0.00000 

-13.22243 

-45.65666 

-14.99417 

-13.22243 

-34.68017 

-25.97066 

-13.22243 

-19.68600 

-29.98834 

-13.22243 

-4.69183 

-25.97066 

-13.22243 

6.28466 

-14.99417 

-13.22243 

10.30234 

» 

-0.00000 

-13.22243 

i  .S .Links .Data  =  [ 

5.91000 

0.00000 

6.89000 

0.00000 

0.00000 

6.89000 

-7.42000 

0.00000 

6.89000 

3; 


if  "isempty (links) ,  D . S . Links . Data  =  links;  end 


HSLJ30XDAT.M 


function  D  =  hsl_boxdat(mass, inertia, links, config, bsize) 


7,  HSL_B0XDAT 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


:  Box  Container  aerodynamic,  mass,  inertial  and  geometric  data 

DAT  =  HSL_B0XDAT( (mass) (.inertia) (.links)  (.config) (,bsize)) 

mass  :  Load  mass  (slugs)  , 

inertia  :  Inertia  matrix  (slug.ft~2) 

links  :  Links  data  matrix  (see  HSL_INIT) 

config  :  Sling  configuration  string  [’single’ I ’multiple’] 

bsize  :  Box  size  [x,y,z]  (ft) 
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l  DAT  :  Data  structure  (see  HSL_INIT) 

7.  R.  A.  Stuckey  01/06/99  (c)  1999,  Defence  Science  and  Technology  Organisation 

7. - 

if  nargin<5,  bsize  =  []  ;  end 

if  nargin<4,  config  =  [] ;  end 

if  nargin<3,  links  =  D;  end 

if  nargin<2,  inertia  =  □;  end 

if  nargin<l,  mass  =  □;  end 

7,  Set  some  defaults 

if  isempty (mass)  ,  mass  =  1000/32.174;  end 

if  isempty (inertia) ,  inertia  =  diag([  0.33  1.20  1.20  ] *mass*32 . 174) ;  end 
if  isempty (bsize) ,  bsize  =  [  6.0  6.0  6.0  ];  end 
if  isempty (config)  ,  config  =  'multiple’;  end 

D  =  struct  ( 'Name' ,  []  ,  'V' ,  [] ,  ’  A’ ,  □  , 'm' ,  []  ,  ’  J’ ,  []  ,  'S’ ,  []  ) ; 

D.Name  =  'BOX’ ; 

D.V  =  [  0.0  ];  7.  Trim  airspeed  (KTAS) 

7,  Construct  the  aerodynamic  coefficient  matrix 

D.A  =  zeros(l, 10*6,1) ; 

7.  v  =  0.0  KTAS 
7. 

A  =  [ 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0. 00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000  0.00000e+000 

]’;  D.  A(1 , :  ,  1)  =  A(:)’; 

D.m  =  mass;  7.  Mass  (slugs) 

D.J  =  inertia;  7.  Inertia  matrix  (slugs. ft"2) 

7.  Construct  the  geometric  struct 

D.S  =  struct  ( 'Patches’ ,  D  , 'Lines’ ,[], 'Links ’,[])  ; 

D.S. Patches. Data  =  [ 

111 
1  1  -1 
-1  1  -1 
-1  1  1 
-1  -1  1 
-1  -1  -1 
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l  -1  -l 

1  -1  l 

] *diag(bsize/2) ; 

D.S. Patches. Coordlndex  =  [ 

0  1  2  3  -1  . . . 

2  3  4  5  -1  . .. 

4  5  6  7  -1  ... 

1  2  5  6  -1  ... 

0  1  6  7  -1  ... 

0  3  4  7  -1 

]; 

if  strcmpi(config, ’single’) 

D.S. Lines. Data  =  [ 

11-1 
0  0-1 
NaN  NaN  NaN 
1  -1  -1 
0  0-1 
NaN  NaN  NaN 
-1  1  -1 
0  0-1 
NaN  NaN  NaN 
-1  -1  -1 
0  0-1 
NaN  NaN  NaN 

] *diag(bsize/2) ; 

D.S. Links. Data  =  [ 

0  0-1 

]*diag(bsize.*[  0.0  0.0  1.0  ]); 

if  "isempty (links) 

if  (length (links)  ==  1) 

D . S. Links .Data(3)  =  links; 
else 

D.S. Links. Data  =  links; 

end 

end 

D. S. Lines. Data([2, 5, 8,11]  , :)  =  ones (4, 1) *D . S .Links .Data; 

elseif  strcmpi(config, ’multiple’) 

D.S. Lines. Data  =  []  ; 

D.S. Links. Data  =  [ 

1  1  -1 
-1  1  -1 
-1  -1  -1 
1  -1  -1 

] *diag(bsize . * [  0.5  0.5  0.5  ]); 

if  "isempty (links) ,  D.S. Links. Data  =  links;  end 

end 
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HSLJLOAD.M 


function  r  =  hsl_load(r ,acj) 


l  HSL.LOAD  :  Calculate  load  positions  from  the  helicopter  positions,  the 
7,  helicopter  and  load  attitudes  and  the  cable  lengths  &  orientation 

I.  r  =  HSL.LOAD  (r  ( ,  ac  j  )  ) 

7. 

7.  r  :  Configuration  position  vector 

7.  =  [  R1*_N  ;  R2*_N  ;  ...  ;  Rn*_N 

7.  al  ;  a2  ■  an 

'/,  acj  :  Cable-angle  matrix  :  Single-cable 

7,  =[[  0  0  0  ]’  ac2  .  . .  acn  ] 


] 

case  only 


7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


The  sub-vectors  are  defined: 

R1*_N  :  helicopter  [x,y,z]  position  vector  in  inertial  axes 
Rj*_N  :  load  [x,y,z]  position  vector  in  inertial  axes 
al  :  helicopter  [phi, theta, psi]  Euler  angular  position  vector 

.aj  :  load  [phi , theta, psi]  Euler  angular  position  vector 

acj  :  cable/cg-line  [phi, theta, psi]  angles  in  cable-axes 


See  HSL.INIT  for  more  information 


7,  R.A.  Stuckey  01/10/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - 


global  HDAT_  LDAT_  CDAT_ 

n  =  length (LDAT_ ) ;  n  =  n+("n); 

if  nargin<2,  acj  =  zeros (n, 3);  end 

'/,  Create  some  transformation  matrices 

T_jN  =  zeros (3,3 ,n) ;  T_cjN  =  NaN*ones(3,3,n) ; 

T_jN(:,:,l)  =  euler(r(n*3+[l:3])) ; 

7,  Determine  sling  configuration 
mj  =  zeros (l,n) ; 

MultipleCables  =  NaN;  SingleCable  =  NaN; 
for  j  =  2:n 

mj(j)  =  length (CDAT_( j) .10) ; 

MultipleCables (j )  =  (mj(j)>2)  k  ("any (diff (sort (CDAT_ (j) . i (1 ,:))))) ; 
SingleCable(j)  =  (mj(j)==l); 

end 

7,  Calculate  the  load  positions 
for  j  =  2:n 

jj  =  ( j-l)*3+[l : 3] ; 


70 


T_jN(:,:,j)  =  euler (r(n*3+j j) ) ;  7.  Load-axes  transformation 

toplinks  =  HDAT_ .S. Links .Data(CDAT_(j) . i (1, : ) , :) ; 
botlinks  =  LDAT_(j)  .S  .Links  .Data(CDAT_(j )  .i  (2, :)  , : )  ; 

if  MultipleCables(j) 

if  any (diff (CDAT_(j) . 1) ) 

error ( ’  Cable  lengths  must  be  equal  ! ’ ) 

end 

'/,  Calculate  load  length  (x)  and  width  (y) 

Rili4_j  =  diff (botlinks ([4,1] ,:))’ ;  b  =  sqrt (sum(Rili4_j .*2)) 
Rili2_j  =  diff (botlinks ( [2,1] ;  c  =  sqrt (sum (Rili2_j.~2)) 

7,  Calculate  distance  from  load  eg  to  helicopter  attachment  point 

zl  =  sqrt (CDAT_ ( j ) . 1 (1) “2- (b/2) “2- (c/2) *2) -botlinks (1 , 3) ; 

7,  Update  the  (load)  position  vector 

r(j j)  =  r(l:3)  ... 

+T_jN( : , : ,1) ’*toplinks (1, :) ’  ... 

+T_jN(:,:,j)’*[  COzl]'; 

elseif  SingleCable( j) 

T_cjN(: , : ,  j)  =  euler (acj  ( : ,  j) )  ;  7.  Cable-axes  transformation 

7.  Update  the  (load)  position  vector 

r(j j)  =  r(l:3)  .  .  . 

+T_jN(: , : ,1) ’*toplinks’  ... 

+T_cjN( : , : , j) ’ * [  0  0  CDAT_(j).l  ]>  ... 

-T_ jN ( : , : , j ) ’ *botlinks ’ ; 

end 

end 


HSLSIM.M 


7.  HSLSIM  :  Helicopter  Slung-Load  dynamic  Simulation 
7. 

7.  HSLSIM 

7. 


7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


Inputs : 

opt_  :  GLOBAL  options  cell  array  [4x1] 

HDAT_  :  GLOBAL  helicopter  struct 
LDAT_  :  GLOBAL  loads  struct  array  [lxn] 

CDAT_  :  GLOBAL  loads  struct  array  [lxn] 
n  :  Number  of  bodies 

t  :  Time  vector  [Nil] 

TD  :  Time  &  control  input  matrix  [Nx5] 

D  :  System  mass/inertia  matrix  [n*6xn*6] 

xO  :  Trim  state  vector  [n*6xl] 
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i  x  :  Initial  state  vector  [n*6xl] 

7. 

’/,  Outputs : 

7. 

7.  U 

7.  R 

7.  R_q 

7.  X 

7.  X_q 

7. 

7.  See  HSL_INIT  for  more  information 

7.  R.A.  Stuckey  01/10/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - 


:  Generalised  velocity  coordinate  matrix  [Nxn*6] 

:  Configuration  position  coordinate  matrix  using  Euler  angles  [Nxn*6] 
:  Configuration  position  coordinate  matrix  using  quaternions  [Nxn*7] 

:  Full  state  variable  matrix  using  Euler  angles  [Nxn*12] 

:  Full  state  variable  matrix  using  quaternions  [Nxn*13] 


global  HDAT_  LDAT_  CDAT_  opt_ 

7.  Check  input  variables  and  set  defaults 

if  length(opt_)<4 

{’combined’ , ’euler’ ,0,1}; 
for  i  =  length(opt_)+l:4 
opt_{i}  =  ans{i}; 

warning ([’  Using  opt_{’ ,int2str(i) , ’}  =  ’,opt_{i}]) 

end 

end 

if  ("strncmpi (opt_{l> , ’longitudinal’ .length (opt_{l})))  ... 

k  ("strncmpi (opt_{l} , ’lateral’ , length (opt_{l>) ) )  ... 

k  ("strncmpi (opt_{l> , ’combined’ , length (opt _{1>) ) ) 

error ( ’  Option,  opt_{l>  =  ’,opt_{l},’  not  recognised  !’) 

end 

opt_2  =  strncmpi (opt_{2>, ’quaternion’ ,length(opt_{2}) ) ;  o2  =  opt_2+3; 
if  ("opt_2)&("strncmpi(opt_{2}, ’euler’ ,length(opt_{2}) ) ) 
error ( ’  Option,  opt_{2>  =  ’,opt_{2},’  not  recognised  !’) 

end 

if  "any(opt_{3}==[0,l]) 

error (’  Option,  opt_{3>  =  ’ , int2str (opt_{3»  ,  ’  not  recognised  ! ’) 

end 

if  opt_{3} 

for  j  =  2:n 

if  isempty(CDAT_(j) .K) 

error(’  Empty  stiffness/damping  matrix,  CDAT_(’ ,int2str(j) , ’) .K  !’) 

end 

end 

end 

if  "any(opt_{4}==[0,l]) 

error ( ’  Option,  opt{4}  =  ’  ,int2str (opt_{4>) , ’  not  recognised  !’) 

end 

if  ("opt_{4})& ("exist (’control’ , ’dir’)) 

error(’  Linear  simulation  requires  CONTROL  SYSTEM  TOOLBOX  !’) 

end 

7.  Determine  sling  configuration 
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mj  =  zeros(l,n) ; 

MultipleCables  =  NaN;  SingleCable  =  NaN; 
for  j  =  2:n 

mj (j)  =  length (CDAT_(j). 10); 

MultipleCables (j)  =  (mj(j)>2)  &  ("any (diff (sort (CDAT_(j) . i (1 ,:))))) ; 

SingleCable ( j )  =  (mj(j)==l); 

end 

7.  Construct  inverse  mass/inertia  matrix 
finD  =  isfinite(diag(D) ) ;  finDn  =  find(finD); 

diag(inf*ones(l,n*6)) ;  ans (finDn, finDn)  =  D (finDn, finDn) ;  D  =  ans; 
zeros(n*6);  ans  (finDn,  finDn)  =  inv(D  (finDn,  finDn))  ;  Di  =  ans; 

7,  Determine  trim  configuration  velocities 

[xdotO.vaO]  =  hsl_fun(0,zeros(n*12,l) ,  []  ,D,Di ,x0 .zeros (n*6, 1)  ,zeros(l,5))  ; 

7.  Set  initial,  perturbed  state  from  trim 

if  "opt_2 
dx  =  x-xO; 
else 

xO_q  =  [x0(l:n*9) ;zeros (n*4,l)] ;  x_q  =  [x(l :n*9) ;zeros(n*4, 1)]  ; 
for  j  =  l:n 

x0_q(n*9+( j-1) *4+[l :4] )  =  e2ql(x0(n*9+( j-1) *3+[l :3] ) ) ; 
x_q(n*9+(j-l)*4+[l  :4]  )  =  e2ql(x(n*9+( j-l)*3+  [1 : 3]  ) ) ; 

end 

dx_q  =  x_q-x0_q; 

end 

7.  Run  the  simulation 

if  opt_{4}==l  7.  Nonlinear  simulation: 

7.  Uses  4/5th  order  R-K  integration  routine  with  fixed  step-size 

options  =  odeset ; 

if  "opt_2  7.  Euler  representation 

[t,dX]  =  ode45f (’hsl_fun’ ,t,dx, options,  D,Di,xO,vaO,TD) ; 

else  7.  Quaternion  representation 

[t,dX_q]  =  ode45f ( ’hsl.fun’ ,t,dx_q, options ,  D,Di,xO_q,vaO,TD) ; 

end 

else  7.  Linear  simulation: 

7.  Uses  linear  time-invariant  time  response  kernel 

if  "opt_2  7.  Euler  representation 

del  =  le-3*ones  (n*12 , 1)  ;  7.7.7.  DELTA  FOR  JACOBIAN  COMPUTATION 

dxO  =  zeros(n*12,l) ; 

7.  Compute  the  state  and  control  jacobian  matrices 
fac  =  []  ;  vec  =  []  ;  del  =  le-3*ones  (n*12 , 1)  ; 

DFDY  =  numjac(’hsl_fun’ ,0,dx0,fty, del, fac, vec,  [],[],[]  ,  D,Di,xO , vaO, zeros (1 ,5) ) ; 
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del  =  le-3*ones(4,l)  ;  l.'/.'i  DELTA 

ydel  =  diag(del);  fdel  =  zeros (n*12,4) ; 
for  j  =  1:4 

fdel( :  ,  j)  =  hsl_fun(0,dx0,  □  ,  D.Di.xO, vaO, [0,ydel( : , j) ’] ) ; 

end 

DFDD  =  (fdel-fty(: .ones (1,4)))  *diag(l  ./del); 

•/.  Convert  linear  system  from  continuous  to  discrete 

sys  =  ss(DFDY,DFDD,eye(n*12) ,zeros(n*12,4)) ; 
sysd  =  c2d(sys,t(2)-t(l) , ’foh’) ; 

l  Implement  any  initial  velocities  and  perform  the  simulation 

rOdot  =  [v0(l:n*3) ;zeros(n*3 , 1)] ; 
for  j  =  l:n 

Wij_j  =  eratesi(x0(n*9+(j-l)*3+[l:3])) ; 

rOdot (n*3+ (j-1) *3+ [1:3] )  =  Wij_j*v0(n*3+(j-l)*3+[l:3])  ; 

end 

dX  =  lsim(sysd,TD(: ,2:5) ,t,dx)+t* [zeros (1 ,n*6) , rOdot ’] ; 

else  '/,  Quaternion  representation 

del  =  le-3*ones  (n*13 , 1) ;  7,7.7.  DELTA  FOR  JACOBIAN  COMPUTATION 

dxO_q  =  zeros(n*13,l) ; 

7.  Compute  the  state  and  control  jacobian  matrices 
fac  =  []  ;  vec  =  []  ;  del  =  le-3*ones (n*12, 1) ; 

DFDY  =  numjac(’hsl_fun’ ,0,dx0_q,fty, del, fac, vec,  [],[],[]  ,  D.Di.xO.q.vaO.zerosCl.S)) 

del  =  le-3*ones(4,l) ;  7.7.7.  DELTA 

ydel  =  diag(del);  fdel  =  zeros (n* 13,4) ; 
for  j  *  1:4 

fdel ( : ,  j)  =  hsl_fun(0,dx0,  []  ,  D,Di,xO,vaO,  [0,ydel(:  ,j) ’])  ; 

end 

DFDD  =  (f del-fty ( : ,ones(l ,4) ) ) *diag(l  ./del); 

l  Convert  linear  system  from  continuous  to  discrete 

sys  =  ss(DFDY,DFDD,eye(n*13) ,zeros(n*13,4)) ; 
sysd  =  c2d(sys,t(2)-t(l) , ’foh’) ; 

7.  Implement  any  initial  velocities  and  perform  the  simulation 

rOdot_q  =  [vO(l :n*3) ;zeros (n*4, 1)]  ; 
for  j  =  l:n 

Wij_j  =  qratesi (x0_q(n*9+( j-1) *4+ [1 :4] ) ) ; 

rOdot _q(n*3+  ( j-1)  *4+  [1:4])  =  Wij_j*vO (n*3+(  j-1) *3+ [1 :3] ) ; 

end 

dX_q  =  lsim(sysd,TD( : ,2:5) ,t ,dx_q)+t*  [zeros (1 ,n*6) , rOdot _q’]  ; 

end 

end 

7.  Create  the  full  state  variable  matrices 
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if  "opt_2 


’/,  Euler  representation 


l  Detect  and  remove  erroneous  (divergent)  observations  from  the  simulation 

isdiv  =  (sign(abs(imag(dX))) lisinf (dX) lisnam(dX))  ; 

if  any (any (isdiv)) 

warning ( ’Divergent  solution  !’) 

N  =  min(find(sum(isdiv,2)))-l;  dX  =  dX(l:N,:); 
else 

N  =  size(dX.l); 

end 

t  =  t(l:N) ;  TD  =  TD(1:N,:); 

l  Add  initial  state  to  perturbed  state  and  calculate  quaternions 

X  =  dX+ones (N , 1) *x0 ’ ;  R  =  X ( : ,n*6+[l :n*6] ) ; 

U  =  X( : , 1 :n*6) ;  R_q  =  [R( 1 :n*3) .zeros (N,n*3)] ; 
for  j  =  l:n 

R_q(:  ,n*3+(j-l)*4+[l  :4]  )  =  e2q(R( :  ,n*3+( j-1) *3+[l :3] ) ) ; 

end 

X_q  =  [U,R_q]  ; 

else  '/,  Quaternion  representation 

'/  Detect  and  remove  erroneous  (divergent)  observations  from  the  simulation 

isdiv  =  (sign(abs(imag(dX_q))) lisinf (dX_q) I isnan (dX_q) ) ; 

if  any (any (isdiv)) 

warning (’Divergent  solution  !’) 

N  =  min(f ind(sum(isdiv,2) ) ) -1 ;  dX_q  =  dX_q(l:N,:); 
else 

N  =  size(dX_q,l) ; 

end 

t  =  t(l:N)  ;  TD  =  TD(1:N, :)  ; 

V,  Add  initial  state  to  perturbed  state  and  calculate  Euler  angles 

X_q  =  dX_q+ones(N,l)*xO_q’ ;  R_q  =  X_q( : ,n*6+ [1 :n*7] ) ; 

U  =  X_q( : , 1 :n*6) ;  R  =  [R_q( : ,l:n*3) ,zeros(N,n*4)] ; 
for  j  =  l:n 

R( :  ,n*3+(j-l) *3+  [1:3] )  =  q2e(R_q( :  ,n*3+( j-1) *4+[l :4] ) ) ; 

end 

X  =  [U,R] ; 

end 

I.  Compute  body-axes  velocities ,  cable  angles  and  internal  cable  forces 

Va  =  zeros(N,n*6) ;  Vadot  =  zeros (N,n*6) ;  TDD  =  zeros(N,5); 
if  n>l 

Acj  =  zeros(N, (n-1) *3) ;  FC  =  zeros(N.n-l) ; 

end 

for  i  =  1:N  '/,  Evaluation  loop  (no  integration) 


l  Call  the  function  with  all  output  arguments 
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if  "opt_2  '/.  Euler  representation 

[xdot , va, vadot , acj  ,fc,d]  =  hsl_fun(t  <i)  ,dX(i ,D ,Di ,x0 ,  vaO ,TD)  ; 

else  7,  Quaternion  representation 

[xdot_q,va, vadot  ,acj  ,fc ,d]  =  hsl_fun(t (i)  ,dX_q(i, :)’,[]  ,D,Di,xO_q,vaO,TD)  ; 

end 

*/,  Insert  the  variables  into  their  respective  observation  matrices 

Va(i , : )  =  va’;  Vadot (i,:)  =  vadot’; 
if  n>l 

squeeze (acj ( : , 1 , 2 : n) ) ;  Acj(i,:)  =ans(:)’; 

FC(i, : )  =  (f c ./diag(D( [2 :n] *3,  [2 :n] *3) ) /32 . 174)  ’ ; 

end 

TDD(i , : )  =  [t  (i) ,d’]  ; 

end 


ODE45F.M 


function  [T,Y,Yd]  =  ode45f (odefile ,T,y0 , options .varargin) 

7,  0DE45F  :  Solve  differential  equations  -  higher  order  method 


7. 

7.  [T,Y,Yd]  =  0DE45F( ’F’ ,T,y0, options ( .varargin) ) 

7. 

•/,  F  :  String  containing  the  name  of  the  ODE  function 

7. 

7.  where  yd  =  F(t,y) 

7. 

7.  T  :  Time  vector  [Nxl] 

7.  yO  :  Initial  conditions  [nxl] 

7. 

7.  Y  :  Solution  matrix  [Nxn] 

7,  Yd  :  Differential  matrix  [Nxn] 

7. 

7.  See  0DE45 ,  ODE 


7.  Modified  from  0DE45.M: 

7. 

7.  C.B.  Moler ,  3-25-87,  8-26-91,  9-08-92. 

7.  Copyright  (c)  1984-94  by  The  MathWorks,  Inc. 

7.  E.A.  Stuckey  10/12/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - 


7.  Initialise  the  integration  parameters  and  matrices 


alf=[  1/4  3/8  12/13  1  1/2]’; 

bet=[  [1  00  0 

[  3  9  0  0 

[  1932  -7200  7296  0 

[  8341  -32832  29440  -845 

[  -6080  41040  -28352  9295 

gam=[  [  902880  0  3953664  3855735 

[  -2090  0  22528  21970 


0 

0 

0 

0 

-5643 

-1371249 

-15048 


0  ]/4 

0  ]/32 

0  ]/2197 

0  ]/4104 

0  ] /20520  ]; 

277020  ] /7618050 
-27360  ] /752400  ]; 
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T  =  T( : ) ; 

N  =  sized, 1);  ny  =  size(yO,l);  Nb  =  size(bet); 

H  =  diff(T);  f  =  zeros(ny,Nb(2)) ;  Y  =  zeros (N,ny);  Yd  =  Y; 

Y(l, :)  =  yO'; 

7,  Set  the  output  function 

if  isempty (options) 
outfun  =  ’ ’ ; 
else 

outfun  =  odeget (options, ’OutputFcn’) ; 

end 

if  "isempty (outfun) 

feval(outfun,T( [1,N] ) ,y0, ’init ’ ) ; 

end 

for  n  =  1:N-1  7.  The  main  loop  . . . 

7,  Compute  the  slopes 

f  ( :  ,  1)  =  feval(odefile,T(n)  ,Y(n, :)’,[]  ,varargin{:})  ; 

Yd(n, :)  =  f(:,l)’; 
for  j  =  l:Nb(l) 

f  ( : ,  j+1)  =  feval(odef  ile,T(n)+H(n)*alf  (j)  ,Y(n, : )  '+H(n)*f  *bet  (j  ,:)’,[]  ,varargin{ : }) ; 

end 

7.  Estimate  the  error 

7. 

7.  delta  =  norm(H(n)*f*gam(2, :)  ’ ,  'inf ')  ; 

7.  Update  the  solution 

Y (n+1 , : )  =  Y(n,:)+H(n)*gam(l,:)*f 
7,  Check  the  status  of  the  output  function 
if  "isempty (outfun) 

status  =  feval(outfun,T(n+l) ,Y(n+l, :)) ; 
if  status,  break,  end 

end 

end 

7.  Evaluate  the  last  time  point  and  truncate  output  matrices  if  necessary 
Yd(n+1 , : )  =  f  eval(odef  ile,  T(n+1)  ,Y(n+l ,:)’,[]  ,varargin{  :>)’ ; 

N  =  n+1;  T  =  T(1:N);  Y  =  Y(1:N,:);  Yd  =  Yd(l:N,:); 

7.  Set  the  output  status 
if  "isempty (outfun) 


77 


DSTO-TR-1257 


feval (outfun,  [],[],  'done’)  ; 

end 


HSL.FUN.M 


function  [xdot ,va, vadot , acj ,Fc ,d]  -  hsl_fun(t ,dx ,FLAG ,D,Di ,x0 , vaO,TD) 


7.  HSL.FUN  :  ODE  function  for  helicopter  slung-load  system 

7. 

7,  [xdot ,va, vadot , acj ,Fc ,d]  =  HSL_FUN(t,dx,FLAG,D,Di ,xO,vaO,TD) 


7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


t  :  Current  time  variable 

dx  :  Current  perturbed  state  vector  [n*6xl] 

FLAG  :  String  with  the  type  of  information  to  return  *N0T  USED* 

D  :  System  mass/inertia  matrix  [n*6xn*6] 

Di  :  Inverse  system  mass/inertia  matrix  [n*6xn*6] 

xO  :  Trim  state  vector  [n*6xl] 

vaO  :  Trim  body-axes  velocity  vector  [n*6xl] 

TD  :  Time  &  control  input  matrix  [Nx5] 

xdot  :  Current  state  rate  vector  [n*6xl] 
va  :  Current  body-axes  velocity  vector  [n*6xl] 
vadot  :  Current  body-axes  acceleration  vector  [n*6xl] 
acj  :  Cable  Jingle  matrix  [3x4xn] 

Fc  :  Internal  cable  force  vector  [n-lxl] 
d  .  :  Current  control  input  vector  [4x1] 

See  HSLSIM  for  details 


7.  R.A.  Stuckey  01/10/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - - - 


global  HDAT_  LDAT_  CDAT_  opt_  tdd_  vadot_  rdot_ 

n  =  length(LDAT_) ;  n  =  n+(~n) ;  nx  =  length(xO); 
opt_2  =  (nx==n*13) ;  o2  =  opt_2+3; 

7.  Early  exit  if  dx  is  in  err 

if  ("isreal(dx) I  any (isinf (dx) lisnan(dx))) 
xdot  =  dx;  return 

end 

7.  Define  some  strings  for  FEVAL  function  evaluation,  below 
if  opt_2 

anglestr  =  'quaternion’;  ratesstr  =  'qratesi'; 
else 

anglestr  =  'euler';  ratesstr  =  'eratesi'; 

end 
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7.  Determine  axes  to  be  used  in  solution  (the  sub-space) 

if  strncmpi (opt_{l} , 'longitudinal' ,length(opt_{l})) 

7.  For  n=2,  ans  =  [1,3,8,4,6,11];  nsys  =  [ans.ans+12] ; 


[1 ,3,n*3+2 ,n*6+ [1 ,3]] ’ tones (1 ,n) +ones (5 , 1) * [0 : n-1] *3 ;  nsys  =  ans(:)’; 
if  opt_2 

[n*9+[l,3]]  ’*ones(l,n)tones(2,l)*[0:n-l]*4; 
else 

[n*9+  [2] ]  ’ tones ( 1 , n)  +one s ( 1 , 1 ) * [0 : n- 1] *3 ; 

end 

nsys  =  [nsys,ans(:)  ’]  ; 

elseif  strncmpi(opt_{l>, ’lateral’ ,length(opt_{l})) 

7.  For  n=2,  ans  =  [2,7,9,5,6,10,12];  nsys  =  [ans,ans+12]  ; 

[2,n*3+[l,3]  ,n*6+[2,3]]  ’*ones(l,n)+ones(5,l)*[0;n-l]*3;  nsys  =  ans(:)’; 
if  opt_2 

[n*9+[l  :4]]  ’ tones (l,n) tones (4, 1) * [0: n-1]  *4; 
else 

[n*9+[l  ,3]]  ’ tones (l,n) tones (2,1)  *  [0:n-l]  *3; 

end 

nsys  =  [nsys,ans(:) ’]  ; 
else 

1  For  n=2,  ans  =  [1:3,7:9,4:6,10:12];  nsys  =  [ans,ans+12] ; 

[1 :3,n*3t [1:3]  ,n*6+[l  :3]  ]  ’ tones (1 ,n) tones (9 , 1) * [0 :n-l] *3 ;  nsys  =  ans(:)’ 
if  opt_2 

[n*9+ [1:4]] ’tones (1 ,n) tones (4 , 1) * [0 : n-1] *4 ; 
else 

[n*9+[l:3]]  ’ tones (l,n) tones (3 , 1)  *  [0:n-l]  *3; 

end 

nsys  =  [nsys ,  cins  ( : )  ’]  ; 

end 

x  =  xOtdx;  7.  Current  state  =  trim  t  perturbation 

7.  Set  all  other  variables,  outside  the  solution  space,  to  zero 

zeros (nx,l);  ans (nsys)  =  x(nsys);  x  =  ans; 
uO  =  x0(l:n+6);  rO  =  x0(n*6+[l:n*(o2+3)]) ; 
u  =  x(l:n*6);  r  =  x(n*6t[l;n*(o2+3)]) ; 

7.  Create  some  transformation  matrices 

T_jN  =  zeros(3,3,n) ;  T.cjN  =  NaN*ones(3,3,4,n) ; 

T_jN(:,:,l)  =  feval(anglestr,r(n*3t[l:o2])) ; 

7,  Determine  sling  configuration 

mj  =  zeros (l,n); 

MultipleCables  =  NaN;  SingleCable  =  NaN; 
for  j  =  2:n 

mj(j)  =  length(CDAT_(j)  .10)  ; 

MultipleCables (j)  =  (mj(j)>2)  &  ('any (diff (sort (CDAT_ (j) . i(l ,:)))))  ; 
SingleCable (j)  =  (mj(j)==l); 

end 

kcj_N  =  NaN*ones (3,4,n) ;  acj  =  NaN*ones(3,4,n) ;  wcj_cj  =  NaM*ones(3,4,n) ; 
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Rlsjs_N  =  NaN*ones (3,n) ; 

'/,  Calculate  the  cable  properties  for  each  load 


for  j  =  2:n 


T_jN(:,:,j)  =  f eval (anglestr ,r (n*3+( j-1) *o2+ [1 : o2] ) ) ;  7.  Load-axes 
Rls js_N( : , j)  =  -r (1 :3)+r ( ( j-1) *3+ [1 :3] ) ; 

toplinks  =  HDAT_ . S .Links .  Data(CDAT_(j) . i (1 ,:),:) ; 
botlinks  =  LDAT_(j) .S. Links. Data(CDAT_(j) .i(2, ; 

if  MultipleCables (j) 

Rajjs_N  =  Rlsjs_N(: , j) *ones (1 ,mj ( j) ) -T_jN ( : , : ,1) ’*toplinks’  ; 

CDAT_(j)  .Raji_N  =  Raj  js_N+T_ jN( : ,  : ,  j )  ’ *botlinks ’ ;  7.  Attachment-attachment 

CDAT_ ( j )  .  1  =  sqrt  (sum( (CDAT_ ( j)  .Raji_N)  .”2,1))  ;  7.  Cable  length 

CDAT_(j)  .ldot  =  u( ( j-2)  *3+3+  [1:3])’*  (T_jN ( :  ,  :  ,j)*CDAT_(j)  ,Raji_N)  . /  (CDAT_(j)  .1); 
elseif  SingleCable (j) 

Rajjs.N  =  Rlsjs_N(: ,j)-T_jN(: , : ,1) ’*toplinks’ ; 

CDAT_(j)  .Raji_N  =  Raj js_N+T_jN( : , : ,  j)  ’*botlinks’ ;  7,  Attachment-attachment 

CDAT_(j).l  =  sqrt(sum((CDAT_(j)  .Raji_N)  ."2,1)) ;  7.  Cable  length 

CDAT_ ( j ) . ldot  =  u( (j-2) *3+3+ [3]); 

end 

for  i  =  l:mj(j) 

kcj_N( : ,i, j)  =  CDAT_ ( j) .Raji_N( : ,i)/CDAT_ ( j) . 1 (i) ; 
ac j ( : , i ,  j)  =  k2a(kcj_N(: ,i,j)) ; 

T_cjN(: , : ,i,j)  =  euler (ac j ( : ,i, j) ) ; 

end 

end 

7.  Construct  the  configuration  matrices 
[A , Ai , Adot ,H]  =  hsl_config(u,T_jN,T_cjN,Rlsjs_N,acj) ; 

l  Correct  the  system  mass/inertia  matrix  for  implementation  into  the  equations 

finD  =  isfinite(diag(D)) ;  finDn  =  find(f inD) ; 

D2  =  zeros (n*6,n*6) ;  D2 (finDn, finDn)  =  D (finDn, finDn) ; 
for  j  =  l:n 

j  j  =  (j-1)  *3+  [1:3]  ; 

D2(j j , j j)  =  [  1  -1  -1  ;  -1  1  -1  ;  -1  -1  1  ] . *D2(jj , j j) ; 

end 

7.  Compute  the  gravitational  force 

g  =  32.174;  g_N  =  [  0  0  g  ]’* [ones (1 ,n) .zeros (1 ,n)] ; 

fg  =  zeros (n*6,l) ; 
for  j  =  l:n 

jj  =  (j-1)  *3+  [1:33  ;  fg(jj)  =  D2(  j  j  (3) ,  j  j  (3) )  *  [  00g]’; 

end 

7.  Calculate  the  body-axes  velocities  and  current  control  inputs 
v  =  A*u;  va  =  v; 
for  j  =  l:n 


7.  Cable  unit-vectors 
7.  Cable  angles 
7.  Transformations 
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jj  =  ( j-l)  *3+  [1:3];  va(jj)  =  T_jN( : , : ,  j)*v(j  j) ; 

end 

if  all(size(TD)==[l,5]) 
d  =  TD(2:5)  ’ ; 
else 

tdi  =  (TD(: ,l)==t) ; 
if  any (tdi) 

d  =  TD (find (tdi) , 2:5) ’ ; 
else 

d  =  interplq(TD( :,l),TD(:,2:5),t)’; 

end 

end 

'/,  Compute  the  aerodynamic  force  and  convert  back  to  inertial  axes 
[AA,BA]  =  hsl.faero(va) ; 

fa  =  AA* (va-vaO)  +  BA*d;  fa  =  diag(D2) .*fa; 
for  j  =  l:n 

jj  =  (j-l) *3+  [1:3];  fa(jj)  =  T_jN(: , : ,  j) ’*fa(j j)  ; 

end 

*/.  Calculate  the  total  static  weight  and  moments  on  the  helicopter 
dl  =  D2(l,l) ;  ml  =  0; 


for  j  =  2:n 

jj  =  (j-l) *3+  [1 : 3]  ; 

dl  =  dl+D2(j j (1) , j j (1)) ;  ml  =  ml+D2( j j (1) , j j (1) )* (rO(l)-rO ( j j (1) ) )  ; 

end 

l  Compute  the  thrust  force  (approximation)  and  convert  back  to  inertial  axes 
ft  =  zeros (n*6, 1) ; 

if  0  ’/.  thrust  (resultant)  is  vertical 

ft (3)  =  -dl*32 . 174; 
ft (n*3+2)  =  -ml*32 . 174; 

else  7,  thrust  is  normal  to  body 

ft  (3)  =  -dl*32. 174; 
ft (n*3+2)  =  -ml*32 . 174; 

for  j  =  l:n 

jj  =  (j-l) *3+ [1:3]  ;  ft(jj)  -  T_jN(: , : ,  j) ’*ft(jj)  ; 

end 

end 

'/.  Compute  the  Coriolis  forces  and  then  the  combined  forces 

X  =  zeros (n*6,l) ; 
for  j  =  l:n 

jj  =  (n+j-l)*3+ [1:3]; 

X(jj)  =  skew3(u(jj))*D2(jj  ,jj)*u(jj)  ; 

end 
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fx  =  -X-D2*Adot*u;  fo  =  fa+fg+ft+fx; 
j,  Determine  a  solution  for  the  cable  forces 
fc  =  zeros (n*6, 1) ;  sf  =  fc; 

H  =  H(finDn,  : )  ;  Di  =  Di (finDn.f inDn) ;  7.  The  sub-space  defined  above 

if  opt_{3}  l  Elastic  solution 

FC  =  zeros(n*12,l) ;  hi  =  zeros (1 ,n*12) ; 

for  j  =  2:n 

jj  =  ( j-1)  *3+  [1 :3]  ; 

if  MultipleCables (j) 

fC  =  (CDAT_(j) . K ( 1 , :)) . *(l-(CDAT_(j) .10) ./ (CDAT_(j) .1))  ... 

+(CDAT_( j) .K(2, :) ) . * (CDAT_( j) . ldot) . / (CDAT_(j) .1) ; 

FC(jj)  =  sum(ones  (3 , 1)  *max (0,f C) .  *  (CDAT_ ( j)  .Raji.N)  ,2)  ;  hi(jj) 
elseif  SingleCable( j ) 

fC  =  (CDAT_ ( j) .K(l) ) . * (l-(CDAT_(j ) . 10) . / (CDAT_ ( j) . 1) )  ... 

+(CDAT_( j) .K(2) ) . *(CDAT_(j) .ldot) ./ (CDAT_(j) . 1) ; 

FC(  j  j  (1) )  =  max(0,fC)  .*(CDAT_(j)  .1) ;  hi(jj(l))  =  1; 

end 

end 

FC  =  FC (find (hi)) ; 
else  7.  Inelastic  solution 

FC  =  -(H' *Di*H)\(H’*Di*fo(finDn)) ; 

end 

fc(finDn)  =  H*FC;  7,  The  cable  forces 

Fc  =  zeros (n-1 , 1) ; 
for  j  =  2:n 

jj  =  ( j  —  1 )  *3+  [1 :3]  ; 

Fc(j-l)  =  sqrt  (f  c  ( j  j)  ’  *f  c  ( j  j) ) ;  7.  The  total  cable  force  for  each  load 

end 

7.  Compute  the  specific  force 
sf(finDn)  =  Di*  (f o  (f  inDn)+f c  (f  inDn)  )  ; 

7,  Calculate  the  accelerations  and  rates  required  by  the  integration 
udot  =  Ai*sf;  rdot  =  zeros (n* (o2+3) , 1) ;  rdot(l:n*3)  =  v(l:n*3); 
for  j  =  l:n 

Wij_j  =  feval(ratesstr,r(n*3+(j-l)*o2+[l:o2] )) ; 

rdot (n*3+( j-1) *o2+ [1 : o2] )  =  Wij_ j*v(n*3+( j-1) *3+ [1 : 3]  )  ; 

end 

7.  Update  the  body-axes  accelerations  for  the  next  iteration 
vdot  =  Adot*u+A*udot ;  vadot  =  vdot; 
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for  j  =  l:n 

jj  =  ( j-1) *3+  [1:3];  vadot(jj)  =  T_jN(: , :  ,j)*vdot(jj) ; 

end 

7,  Combine  the  velocities  and  positions  into  the  state  rate  vector 
xdot  =  [udot ; rdot] ; 

7.  Augment  the  state  rate  with  the  null  axes  (non-solution  space) 
zeros (nx,l);  ans(nsys)  =  xdot(nsys);  xdot  =  ans; 


HSLJFAERO.M 

function  [A,B]  =  hsl_faero(va) 

7.  HSL_FAER0  :  Function  to  compute  aerodynamic  state-space  matrices 

7. 

7.  [A,B]  =  HSL_FAERO(va) 

7. 

7.  va  :  Body-axes  velocity  vector  [n*6xl] 

7. 

7  A  :  State  matrix  [n*6xn*6] 

7,  B  :  Control  matrix  [n*6x4] 

7. 

I  See  HSL_FUN  for  more  information 

7.  R.A.  Stuckey  01/10/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - 

global  HDAT_  LDAT_  opt_ 
n  =  length (va)/6; 

7.  Determine  aerodynamic  derivatives  to  be  used  in  solution 

if  strncmpi(opt_{l>, ’longitudinal’ ,length(opt_{l})) 

[10  10  10  ]’*[  1010101001];  1  =  ans ( : ) ’  ; 
elseif  strncmpi(opt_{l>,  ’lateral’  ,length(opt_-(l})) 

[010101  ]’*[  0101010111];  1  =  ans  ( : )  ’ ; 
else 

I  =  ones(l,60); 

end 

7.  Aerodynamic  coefficient  matrices  for  the  helicopter 

Vj  =  sqrt  (va(l :  3)  ’  *va(l  :3) )  *0. 5925;  nv  =  length (HDAT_ -V)  ;  7.  Velocity  (knots) 
if  Vj<HDAT_.V(l) 

Vj  =  HDAT_ . V(l) ; 
elseif  Vj>HDAT_.V(nv) 

Vj  =  HDAT_ . V (nv) ; 

end 

if  nv>l 

7.  Interpolate  between  state/control  matrices  based  on  velocity 
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if  1  7.7.7.  SAS  off 

Aj  =  interplq(HDAT_.V,HDAT_.A(: , : ,1) ,Vj) .*1; 
else  7.7.7.  SAS  on 

Aj  =  interplq(HDAT_ . V,HDAT_. A( : , : ,2) , Vj) . *1 ; 

end 

else 

Aj  =  HDAT_ . A. *1 ; 

end 

7.  Integrate  derivatives  into  the  system  state-space  matrices 

a  =  zeros(6,6);  A  =  zeros(n*6,n*6) ; 
b  =  zeros(6,4);  B  =  zeros (n*6 ,4) ; 

jj  =  [l:3,n*3+[l:3]]  ; 

a( : )  =  Aj (1:36) ;  A(jj.jj)  =  a; 
b( : )  =  Aj  (37 : 60)  ;  B(jj,:)  =  b; 

7.  Aerodynamic  coefficient  matrices  for  each  load 

for  j  =  2:n 

jj  =  ( j-1) *3+ [1:3] ; 

Vj  =  sqrt (va( j j) ’ *va( j j) ) *0. 5925 ;  nv  =  length(LDAT_(j) . V) ; 
if  Vj<LDAT_(j) .V(l) 

Vj  =  LDAT_ ( j ) . V(l) ; 
elseif  Vj>LDAT_(j) -V(nv) 

Vj  =  LDAT_(j)  .VCnv) ; 

end 

7,  Interpolate  between  state/control  matrices  based  on  velocity 
if  nv>l 

Aj  =  interplq(LDAT_(j) .V,LDAT_(j) .A,Vj) .*1; 
else 

Aj  =  LDAT_(j) . A.*I; 

end 

7.  Integrate  derivatives  into  the  system  state-space  matrices 

jj  =  [(j-1)  *3+  [1:3]  ,n*3+(  j-1)  *3+  [1 :3]  ]  ; 

a( : )  =  Aj (1:36) ;  A(jj,jj)  =  a; 
b(  :  )  =  Aj  (37 :  60)  ;  B(jj.:)  =  b; 

end 

if  0  7.7.7.  CONTROL  INPUTS  ARE  IN  0.1  IN  UNITS 
B  =  B+0.1; 

end 


HSL_CONFIG.M 

function  [A,Ai,Adot,H]  =  hsl.conf ig(u,T_jN,T_cjN,Rls js_N ,ac j ) 

7.  HSL_CONFIG  :  Calculate  the  configuration  matrix  and  its  inverse 


84 


DSTO-TR-1257 


7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


[A,Ai,Adot,H]  =  HSL_CONFIG(u,T_jN,T_cjN,Rlsjs_N,acj) 

u  :  Generalised  velocity  vector  [n*6xl] 

T_jN  :  Load-axes  transformation  matrix  [3x3xn] 

T_cjN  :  Cable-axes  transformation  matrix  [3x3xn] 

Rlsjs.N  :  Position  vector  from  helicopter  to  load  eg  [3xn] 
acj  :  Cable  angle  matrix  [3x4xn] 

A  :  Configuration  matrix  [n*6xn*6] 

Ai  :  Inverse  configuration  matrix  [n*6xn*6] 

Adot  :  Configuration  rate  matrix  [n*6xn*6] 

H  :  Basis  matrix  [n*6xhh] 

See  HSL_FUN  for  more  information 


7.  Also  see  reference  [1]  for  details  of  the  configuration  matrices. 

7. 

7.  1.  Stuckey,  R.A 

7,  "Mathematical  Modelling  of  Helicopter  Systems" 

7.  DSTO-TR-OOO,  Defence  Science  and  Technology  Organisation,  Sep,  1999 
7.  R.A.  Stuckey  28/07/98  (c)  1998,  Defence  Science  and  Technology  Organisation 

7. - 


global  HDAT_  LDAT_  CDAT_  opt. 
n  =  length (LDAT.) ;  n  =  n+(~n); 

7,  Determine  sling  configuration 
mj  =  zeros (1 ,n) ; 

MultipleCables  =  NaN;  SingleCable  =  NaN; 
for  j  =  2:n 

mj(j)  =  length(CDAT_(j)  .10)  ; 

MultipleCables (j)  =  (mj(j)>2)  &  (*any(diff (sort(CDAT_ (j) .i(l, :))))) ; 
SingleCable  (j)  =  (mj(j)==l); 

end 

7.  Construct  the  configuration  matrices 

A  =  eye(n*6,n*6) ;  Ai  =  A;  Adot  =  zeros (n*6,n*6) ;  H  =  zeros(n*6,n*12) ; 


for  j  =  2:n 

jj  =  ( j-1)  *3+  [1 : 3]  ; 

toplinks  =  HDAT_. S. Links. Data(CDAT_(j) .i(l, ; 
botlinks  =  LDAT. (j)  . S .Links  .DataCCDAT. (j ). i  (2 ; 

if  MultipleCables (j) 

Rjsaj_j  =  T_jN( : , : , j) * (-Rls js_N( : , j)+T_jN(: , : ,1) ' *toplinks (1 , :) ’) ; 

A(j  j  ,1:3)  =  eye(3) ; 

A(jj,  jj)  =  T_jN(: , : , j) ’ ; 

A(  j  j  ,n*3+  [1:3])  =  -T_  jN  C :  ,  : ,  1)  ’  *skew3  (toplinks  (1 , :  ) )  ;  7.  A_j,n+1 
A(j  j  ,n*3+j j)  =  T_jN( :  , : ,  j)  ’ *skew3(Rjsaj_j)  ;  7.  A.j  ,n+j 
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if  nargout>l 

Ai(j j ,1:3)  =  -T_jN(: , : , j)  ; 

Ai(jj , j j)  =  T_jN ( : , : , j) ; 

Ai  ( j  j  ,  n*3+  [1 : 3]  )  =  -T_jN( : , :  ,  j )  *A(  j  j  ,n*3+  [1 : 3]  ) ;  7.  B_j,n+1 

Ai(j  j  ,n*3+j j)  =  -T_jN(: ,  :  ,  j)*A(jj  ,n*3+jj)  ;  7.  B_j,n+j 

end 

if  nargout>2 

Adot(jj.jj)  =  T_jN( : , : ,  j)  ’*skew3(u(n*3+j j ) )  ; 

Adot  (j  j  ,n*3+  [1 : 3]  )  =  ... 

-T_jN(: , :  ,1)  ’*skew3(u(n*3+[l  :3] ) )  *skew3(toplinks(l , : ) )  ;  7.  Adot_j,n+l 

Adot (j j ,n*3+j j)  =  ... 

T_jN(: , : ,  j) ’*(  skew3(u(n*3+j  j))*skew3(Rjsaj_j)-skew3(u(jj))  );  7.  Adot_j,n+j 

end 

if  nargout>3 

H(l:3, j j)  =  eye (3) ; 

H(jj,jj)  =  -eye  (3); 

H(n*3+[1:3]  ,  jj)  =  skew3(toplinks (1 , : ) ) *T_jN( : , : , 1) ;  7.  Zlj_l 
H(n*3+j j  ,  j  j)  =  -skeu3(Rjsaj_j)*T_jN(: , :  ,j);  7.  Zjj_j 

end 

elseif  SingleCable(j) 

A( j j , 1 : 3)  =  eye (3) ; 

A(jj.jj)  =  T_cjN( : , : , 1 ,  j)  ’ ; 

A( j j ,n*3+ [1:3])  =  -T_jN( : , : , 1) ’ *skew3(toplinks) ;  7.  A_j,n+1 

A(jj,n*3+jj)  =  T_jN(: ,  :  ,j)  ’ *skew3(botlinks)  ;  7.  A_j  ,n+j 

if  nargout>l 

Ai( j j ,1:3)  =  -T_cjN ( : , : , 1 , j)  ; 

Ai  ( j  j  ,  j j)  =  T_cjN(: , : ,1, j)  ; 

Ai  (j  j  ,n*3+[l  :3]  )  =  -T_cjN(: , :  ,l,j)*A(jj,n*3+[l:3])  ;  7.  B_j,n+1 

Ai(j  j  ,n*3+j  j)  =  -T_cjN(:,:,l,j)*A(jj,n*3+jj);  7.  B_j,n+j 

end 

if  nargout>2 

acjdot  =  hsl_config_ac(acj(l,l,j) ,CDAT_(j) .l,u(jj))  ; 

Wcj_cj  =  erates (ac j ( : , 1 , j ) ) ;  wcj_cj  =  Wcj_cj*acjdot ; 

Adot (j j , j j)  =  T_cjN(:,:,l,j)’ »skew3(wcj_c j) ; 

Adot(jj  ,n*3+[l:3]  )  =  ... 

-T_jN(: , : ,  1)  1  *skew3  (u(n*3+  [1 :3]  ) )  *skew3(toplinks)  ;  7.  Adot_j  ,n+l 

Adot (j j ,n*3+j j)  =  ... 

T_jN(:  , :  ,j)  ’*(  skew3  (u(n*3+j  j ) )  *skew3(botlinks)  );  7.  Adot_j,n+j 

end 

if  nargout>3 

H(l:3, j j)  =  eye (3) ; 

H(jj.jj)  =  -eye (3)  ; 

H(n*3+[1:3]  ,  jj)  =  skew3(toplinks)*T_jN(:  , :  ,1)  ;  7.  Zlj_l 

H(n*3+j  j  ,  j  j)  =  -skew3(botlinks)*T_jN(:  , :  ,j)  ;  7.  Zjj_j 

end 

end 

if  0  7.7.7.  CHECK  EXPLICIT  INVERSE 

fprintf (’\n  CHECKING  INVERSE  ...’) 
inverr  =  mai(max(abs(inv(A)-Ai)))  ; 
if  (inverr>le-6) 

error (sprintf(’  Inverse  not  accurate:  ERR  =  7.g’  .inverr)) 

end 

end 

end 
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'/.  Construct  the  basis  matrix 
hi  =  zeros(l,n*12) ; 

kcj_N  =  zeros(3,n);  kcl_N  =  kcj_N;  kc2_N  =  kcj_N; 

for  j  =  2:n 

jj  =  ( j-1) *3+  [1 : 3]  ; 

if  MultipleCables  ( j )  '/,  Use  entire  force-space 

hi(jj)  =  1; 

elseif  SingleCable(j)  ’/.  Use  only  cable  z-axis 

kcj_N(:,j)  =  T_cjN(3, : ,1, j) ’ ; 

H( : ,  j  j  (D)  -  H(:,jj)*kcj_N(:,j);  hi(jj(D)  =  1; 

end 

end 

H  =  H(: .find(hi)); 


7. - 

function  acdot  =  hsl_conf ig_ac (ac ,1 ,vc) 
l  hsl_config_ac  :  Calculate  cable  angular  rates 

7. 

7.  acdot  =  HSL_CONFIG_AC  Cac ,  1 , vc) 

7. 

%  ac  :  Cable  angle  vector  [lxn] 

’/  1  :  Cable  length  vector  [lxn] 

\ l  vc  :  Cable  velocity  vector  [lxn] 

7. 

7,  acdot  :  Cable  angle-rate  matrix  [3xn] 

7. 

7.  See  HSL.CONFIG 

m  =  length (1) ; 

acdot  =  [  vc(l) ,/l./cos(ac(l, :))  ;  -vc(2)./l  ;  zerosCl.m)  ]; 


EULER.M 

function  T  =  euler(e) 

7.  EULER  :  Euler-angle  transformation  matrix 

7. 

7.  T  =  EULER(e) 

7. 

7.  e  :  vector  of  Euler  angles  [1x3] 

7.  =  [  phi  theta  psi  ] 

7. 

7.  T  :  transformation  matrix  [3x3] 

7. 

7,  where  the  transf ormation  from  inertial  (N)  to  body  (b)  axes  follows: 
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t 

•/.  F_b  =  T_bN  *  F_N  ;  T_bN  =  euler(e) 

7. 

•/.  See  HSL_FUN 

'/.  R.A.  Stuckey  01/10/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - 


ce  =  cos(e);  se  =  sin(e)  ; 


T  =  [  ce(2)*ce(3) 

se(l)*se(2)*ce(3)-se(3)*ce(l) 
ce  (3)  *ce  (1)  *se  (2)+se  (3)  *se  (1) 


ce(2)*se(3) 

se(3)*se(2)*se(l)+ce(3)*ce(l) 

se(3)*ce(l)*se(2)-ce(3)*se(l) 


-se(2) 
se(l)*ce (2) 
ce(l)*ce(2)  ]; 


E2Q.M 

function  Q  =  e2q(E) 

f  E2Q  :  Convert  angular  positions  from  Euler  to  quaternion  representation 

7. 

7.  Q  =  E2Q(E) 

7. 

7,  E  :  Euler  matrix  [mx3]  I  [3xn] 

7,  =  [  Phi  Theta  Psi  ] 

7. 

7,  Q  :  quaternion  matrix  [mx4]  I  [4xn] 

7.  =  [  E0  El  E2  E3  ] 

7. 

7.  See  also  Q2E 

7,  R.A.  Stuckey  07/05/98  (c)  1998,  Defence  Science  and  Technology  Organisation 

X - 

Ne  =  size(E);  trans  =  (Ne(2)"=3); 

if  trans,  E  =  E’ ;  Ne  =  Ne([2,l]);  end 

SE  =  sin (E/2) ;  CE  =  cos(E/2);  Q  =  zeros (Ne (1) ,4) ; 

Q(:,l)  =  SE( : ,1) . *SE( : ,2) . *SE( : ,3)  +  CE( : , 1) . *CE( : ,2) . *CE( : ,3) ; 

Q(:,2)  =  -CE( : , 1) . *SE( : ,2) . *SE( : ,3)  +  SE( : , 1) . *CE( : ,2) . *CE( : ,3) ; 

Q(:,3)  =  SE( : , 1) . *CE( : , 2) . *SE( : ,3)  +  CE ( : , 1) . *SE( : ,2) . *CE ( : ,3) ; 

Q ( : ,4)  =  CE( : , 1) . *CE( : , 2) . *SE( : ,3)  -  SE( : , 1) . *SE( : ,2) . *CE( : ,3) ; 

if  trans ,  Q  =  Q ’ ;  end 


ERATES.M 

function  W  =  erates(e) 

7.  ERATES  :  Euler  transform  matrix  for  angular  velocities 

7. 

7.  W  =  ERATES  (e) 
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7. 

7.  e  :  vector  of  Euler  angles  [1x3] 

’/,  =  [  phi  theta  psi  ] 

7. 

7.  W  :  transf ormation  matrix  [3x3] 

7. 

7.  where 

7.  ... 

7,  [pqr]’=W*[  phi  theta  psi  ]  ’ 

7. 

7.  See  also  ERATESI 

7.  R.A.  Stuckey  07/05/98  (c)  1998,  Defence  Science  and  Technology  Organisation 

7. - 

ce  =  cos(e);  se  =  sin(e)  ; 

W  =  [  1  0  -se (2) 

0  ce(l)  se(l)*ce(2) 

0  -se(l)  ce(l)*ce(2)  ]; 


ERATESI.M 


function  Wi  =  eratesi(e) 

7.  ERATESI  :  Inverse  Euler  transform  matrix  for  angular  velocities 

7. 

7.  Wi  =  ERATESI  (e) 

7. 

7,  e  :  vector  of  Euler  angles  [1x3] 

7,  =  [  phi  theta  psi  ] 

7. 

7,  Wi  :  inverse  transformation  matrix  [3x3] 

7. 

7.  where 

7.  ... 

7.  [  phi  theta  psi  ]’  =  Wi  *  [  p  q  r  ]  ’ 

7. 

7.  See  also  ERATES 

7.  R.A.  Stuckey  07/05/98  (c)  1998,  Defence  Science  and  Technology  Organisation 

7. - 


ce  =  cos(e);  se  =  sin(e) ;  te  =  tan(e) ; 


Wi  = 


[  1 
0 
0 


se (1) *te (2) 
ce(l) 

se(l)/ce(2) 


ce(l)*te(2) 

-se(l) 

ce(l)/ce(2)  ]; 


K2A.M 

function  acj  =  k2a(kcj) 
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'/.  K2A  :  Convert  k-unit  vectors  to  cable-oriented  Euler  angles 

7. 

7.  a  =  K2A(k) 

7. 

i  k  :  unit  vectors  [nx3] 

7. 

7.  a  :  Euler  angle  vector  [nx3] 

7. 

7,  See  HSL_FUN 

7.  R.A.  Stuckey  20/04/99  (c)  1999,  Defence  Science  and  Technology  Organisation 

7. - 

nk  =  size(kcj); 

if  -any(nk==3) ,  error (’  size(kcj ,2)"=3! ’) ,  end 
if  nk(2)"=3,  kcj  =  kcj’;  end 
acj  =  0*kcj; 

acj (: ,1)  =  asin(-kcj (: ,2)) ; 

acj(:,2)  =  asin(  kcj  ( : ,  1) . /cos  (acj  ( : ,  1)  )  ); 

acj  (:  ,3)  =  0; 

if  nk(2)"=3,  acj  =  acj';  end 


SKEW3.M 

function  S  =  skew3(v) 

7.  SKEW3  :  The  general  skew-symmetric  matrix,  S(x,y,z) 

7. 

7.  S  =  skew3(v) 

7. 

7,  v  :  Generic  [x,y,z]  vector 

7. 

7,  S  :  Skew  matrix  [3x3] 

7, 

7.  See  HSL.CONFIG,  HSLJFUN 

7.  R.A.  Stuckey  20/04/99  (c)  1999,  Defence  Science  and  Technology  Organisation 

7. - 

S  =  [  0  -v(3)  v(2) 

v (3)  0  -v(l) 

-v(2)  v(l)  0  ]; 


HSLPLOT.M 

7.  HSLPLOT  :  Plot  output  from  HSLSIM 

7. 

7.  HSLPLOT 
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7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


Inputs: 

opt.  :  GLOBAL  options  cell  array  [4x1] 


Va 

Vadot 

FC 

E 

Acj 

TDD 


Body-axes  velocity  matrix  [Nxn*6] 

Body-axes  acceleration  matrix  [Nxn*6] 

Cable  force  matrix  [Nx(n-l)] 

Configuration  position  coordinate  matrix  [Nxn*6] 
Cable  angle  matrix  [Nx(n-1)*3] 

Updated  time  ft  control  input  matrix  [Nx5] 


See  HSL.INIT,  HSLSIM 


7.  R.A.  Stuckey  01/10/97  (c)  1997,  Defence  Science  and  Technology  Organisation 

7. - 


yv  =  {’Rotation’ ,0, ’Horizontal Alignment’ , ’right’}; 

7,  Determine  axes  to  be  plot 

if  strncmpi (opt_{l}, ’longitudinal’ , length (opt_{l>) ) 


Plots  =  {  ’u’  ’t’  ;  ’w’  ’tc’  ; 
if  "any(any (TDD(xn, [2,5]))) 

;  ’q’ 

’db’  ; 

;  ’fc’ 

’dc’  > 

Plots  =  {  ’u’  ’t’  ;  ’w’  ’tc’  ; 
end 

;  ’q’ 

’ud’  ; 

;  ’fc’ 

’wd’  } 

elseif  strncmpi (opt_{l>, ’lateral’ ,length(opt_{l})) 


Plots  =  {  ’v’ 

’h’  ;  ’w’ 

’t’  ; 

’r ’  ’s’ 

;  ’fc’ 

’da’  ;  ” 

’dr’  ; 

’  ’dc 

if  "any(any (TDD(xn, [3:5]))) 

Plots  =  {  ’v’ 

’h’  ;  ’w’ 

’t’  ; 

>r’  >s> 

H, 

O 

’vd’  ;  ” 

’wd’  >; 

end 

Plots  =  {  ’v’ 

’h’  ;  ’w’ 

’t’  ; 

’r ’  ’s’ 

;  ’fc’ 

’vd’  >; 

elseif  strncmpi(opt_{l}, ’combined’ , length (opt_{l>)) 


Plots  =  {  ’u’  ’p’  ;  ’v’  ’q’  ; 
if  "any (any (TDD (xn,  [2:5]))) 

’  w  ’ 

’r’  ; 

’h’  ’da’  ; 

’he’ 

’db’ 

;  ’f 

’dr’  ; 

;  ’tc 

Plots  =  {  ’u’  ’h’  ;  ’v’  ’he’ 
end 

;  ’w’ 

;  ’p’  ’tc’ 

:  ’q’ 

’ud’ 

.  ,r> 

’vd’  ; 

;  ’fc 

end 

sp  =  {size (Plots,  1)  ,size(Plots,2))-; 

Plots (sp{l}+l, : )  =  {  ”  ”  >;  sp{l>  =  sp{l>+l; 

if  (sp{l}>8) ,  warning ( ’  Too  many  subplots  !  Figure  size  will  be  reduced.’),  end 
ha  =  zeros (sp{ :}) ;  hL  =  zeros(n,sp{l}*sp{2}) ; 
pos  =  get (0 , 'Def aultFigurePosition’ ) ; 

pos(2)  =  pos(2)-pos(4)*(sp{l}-4)/4;  pos(4)  =  pos(4)*sp{l}/4;  pos(3)  =  pos(3); 
h  =  figure( ’Position’ , pos) ; 

set (h, ’Def aultAxesPosition’ , [0.26  0.19  0.70  0.73]) 


’dc’  }; 
’wd’  }; 
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for  i  =  l:sp{l> 
for  j  =  1 : sp{2} 

if  "isempty(Plots{i, j}) 
k  =  (i-1) *sp{2}+j ; 
ha(i,j)  =  subplot (sp{ : },k) ; 

if  (i<sp{l})& ("is empty (Plots{i+l, j})) 

set (ha(i , j) , ’XTickLabel’ ,  □  ) ; 
else 

xlabeK’t,  sec’) 

end 

end 

7,  Plot  each  set  of  variables 

nn  =  [0:n-l]  ; 
switch  Plots{i,j} 
case  ’u’ 

hL(:,k)  =  plot(t(xn)  ,Va(xn,nn*3+i) , ’-k’) ;  ylabel({’u’ , ’ft/s’},yv{:}) 
case  ’v’ 

hL(:,k)  =  plot (t (xn) , Va(xn,nn*3+2) , ’-k’ ) ;  ylabel ({’ v’ , ’ft/s ’} ,yv{ : }) 
case  ’w’ 

hL(:,k)  =  plot(t(xn) ,Va(xn,nn*3+3) , ’-k’) ;  ylabel({’w’ , ’ft/s ’} ,yv{ :}) 
case  ’p’ 

hL(:,k)  =  plot(t(xn) ,Va(xn,n*3+nn*3+l) *180/pi , ’-k’) ;  ylabel(-C’p’ , ’deg/s ’},yv{:}) 
case  ’q’ 

hL(:,k)  =  plot(t(xn) ,Va(xn,n*3+nn*3+2) *180/pi ,  ’-k’)  ;  ylabel ({ ’q’ , ’deg/s ’ },yv{ : }) 
case  ’r’ 

hL(:,k)  =  plot(t(xn) ,Va(xn,n*3+nn*3+3) *180/pi , ’-k’) ;  ylabel ({ ’r’ , ’deg/s ’> ,yv{ :}) 
case  ’fc’ 

hL(2:n,k)  =  plot(t(xn)  ,FC(xn, : ) ,  ’  — k ’ )  ;  ylabel({’fc’  ,  ’/g’},yv{:» 
case  ’h’ 

hL(:,k)  =  plot  (t  (xn)  ,R(xn,n*3+nn*3+l>*  180/pi ,  ’  -k  ’ ) ;  ylabel  ({’\phi’ ,  ’deg’ }  ,yv-[ : » 
case  ’t’ 

hL(:,k)  =  plot (t (xn) ,R(xn,n*3+nn*3+2) * 180/pi ,  ’  -k ’ ) ;  ylabel ({’ \theta’ , ’deg’ > ,yv{ : }) 
case  ’s’ 

hL(:,k)  =  plot(t(xn) ,R(xn,n*3+im*3+3)*180/pi, ’-k’) ;  ylabel ({’\psi’ , ’deg’},yv{ :}) 
case  ’he’ 

hL(2:n,k)  =  plot (t (xn) , Acj (xn, [0 :n-2] *3+l)*180/pi, ’-k’ ) ; 

ylabel ({’\phi_c’ , ’deg’},yv{:» 
case  ’tc’ 

hL(2:n,k)  =  plot (t (xn) , Acj (xn,  [0 : n-2] *3+2) *180/pi, ’-k’ ) ; 

ylabel ({’\theta_c’ , ’deg’} ,yv{ : }) 
case  ’sc’ 

hL(2:n,k)  =  plot(t(xn) ,Acj (xn, [0:n-2]*3+3)*180/pi, ’-k’) ; 

ylabel ({’\psi_c’ , ’deg’ } ,yv{ : }) 
case  ’db’ 

hL(l,k)  =  plot(t(xn) ,TDD(xn,2)) ;  ylabel({ ’ \delta_b’ , ’in’} ,yv{ : }) 
case  ’da’ 

hL(l,k)  =  plot(t(xn) ,TDD(xn,3)) ;  ylabel ({’\delta_a’ , ’in’},yv{:}) 
case  ’dr’ 

hL(l,k)  =  plot(t(xn) ,TDD(xn,4)) ;  ylabel ({’\delta_r’ , ’in’},yv{:}) 
case  ’dc’ 

hL(l,k)  =  plot(t(xn) ,TDD(xn,5)) ;  ylabel ({’ \delta_c ’ , ’in’},yv{:}) 
case  ’ud’ 

hL(:,k)  =  plot(t(xn) ,Vadot (xn,nn*3+l) , ’-k’) ; 

ylabel ({ ’ \f ontsize{6}\bullet ’ , ’ \f ontsize{10}u’ , 'f t / s"2 ’} ,yv{: }) 
case  ’vd’ 
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hL(:,k)  =  plot(t(xn) ,Vadot(xn,nn*3+2) , ’-k’) ; 

ylabel({’\fontsize{6}\bullet’ , ’\f ontsize{10}v’ , ’ft/s“2’},yv{:}) 
case  ’wd’ 

hL(:,k)  =  plot(t(xn) ,Vadot(xn,nn*3+3) , ’-k’) ; 

ylabel({’\fontsize{6}\bullet’ , ’ \f ontsize{10}w’ , ’ft/s*2’},yv{:}) 
case  ’pd’ 

hL(:,k)  =  plot(t(xn) ,Vadot(xn,n*3+nn*3+l)*180/pi, ’-k’) ; 
ylabel({’\fontsize{6)Abullet’ ,  ’\f  ontsize{10}p’ ,  ’deg/s"2’},yv{:]-) 
case  ’qd’ 

hL(:,k)  =  plot  (t  (xn)  ,Vadot  (xn,n*3+nn*3+2)*180/pi ,  ’-k’ ) ; 
ylabel({’\fontsize{6}\bullet’ ,  ’\fontsize{10]-q’ ,  'deg/s~2’},yv{:}-) 
case  ’rd’ 

hL(:,k)  =  plot(t(xn) ,Vadot(xn,n*3+nn*3+3)*180/pi, ’-k’) ; 
ylabel({’\fontsize{6}\bullet’ , ’ \f ontsize{10}r ‘ , ,deg/s"2’},yv{:}) 

end 

end 

end 

SP  =  diag(  [1 :  sp{l}]  )  *  (ha>0) ;  [ans.ans]  =max(SP);  [spi.spj]  =max(ans); 

set(ha(find(ha))  ,  ’XLim’ ,  [0  round (t(xn(end)))]) 
set (ha(find(ha)) , ’XGrid’ , ’on’ , ’YGrid’ , ’on’) 

’/,  Change  line  properties  for  clarity 

set(hL(l,find(hL(l, :))), ’LineStyle’ ’Linewidth' ,1.0, ’Color’ , [  000]) 

if  n>l,  set(hL(2,find(hL(2,:))), ’LineStyle’,’-’, ’LineVidth’ ,0.5, ’Color’, 0.5*C  1  1  1 
if  n>2,  set(hL(3,  find(hL(3, :))) ,  ’LineStyle’,’  —  ’,  ’LineWidth’,  0.5,  ’Color’,  0.5*[  1  1 
if  n>3,  set (hL(4,find(hL(4, :))), ’LineStyle’ ’LineWidth’ ,0.5, ’Color’ ,0.5*[  111 

’/.  Add  a  legend 

if  n==l 
return 
elseif  n==2 

[hi, ho]  =  legend(hL(l :2, 1) , ’Helicopter’ , ’Slung-load’ ,2) ; 
elseif  n==3 

[hi, ho]  =  legend(hL(l :3, 1) , ’Helicopter’ , ’Aft  Load’ , ’Forward  Load’, 3); 
elseif  n==4 

[hi, ho]  =  legend(hL(l:4,l) , ’Helicopter’ , ’AftLoad’ , ’MidLoad’ , ’FwdLoad’ ,4) ; 

end 

7.  Remove  those  pesky  callbacks  (created  by  LEGEND) 

set (hi, ’ButtonDownFcn’ , ’ ’ ,  ... 

’DeleteFcn’ ,  ”  ,  ... 

’UserData’ ,  [] ) 
set (ho, ’ButtonDownFcn’ ,  ’  ’) 

7,  Fix  up  the  legend  text 

ht  =  f indob j (ho, ’Type’ , ’text’ ) ; 

set(ht,’FontUnits’ , ’points’) ; 

tpl  =  get (ht (1) , ’Position’) ; 
set (ht(l) , 'Units’ , ’pixels’) 


] ) ,  end 
1  ] ) ,  end 
] ) ,  end 
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a5e  =  get(ht(l) , 'Extent’) ; 

set  (ht  (1) , 'Units’ , 'data' , 'Position' ,tpl) 

set  (hi, 'Units’ , ’pixels’) 
lp  =  get (hi, 'Position') ; 
lp(4)  =  1.5*a5e(4); 
set(hl, 'Position' ,lp) 
set  (hi, 'Units’ .’normalized’) 

lp  =  get (hi, 'Position') ; 

alp  =  get(ha(l,l) , 'position')  ;  a2p  =  get (ha(l ,2) , 'position' ) ; 

a3p  =  get(ha(spi-l,spj) , 'position') ;  a4p  =  get(ha(spi,spj) , 'position') ; 

lp(l)  =  (alp(l)+alp(3)+a2p(l) )/2-lp(3)/2; 
lp(2)  =  a4p(2)-2*(a3p(2)-(a4p(2)+a4p(4)))-lp(4); 
set  (hi, 'Position' ,lp) 
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